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SECTION 1 
INTRODUCTION 

1 .1 Simulation Extensions 

Within the context of the overall ACOSS Eleven Program, the objec¬ 
tive of the Simulation Extensions Project has been to identify and quantif 
those generic mission, scenario and sensor parameters (as well as their 
interactions) which drive the performance of staring infrared ( IR) space- 
based surveillance systems. This project has been divided into two broad 
tasks: Simulation Enhancements, and Simulation Applications. 

Work in the Simulation Enhancement Task has been directed toward 
generalizing those models in the Draper Integrated Simulations (DIS) which 
are too limited in scope or perhaps restricted to one particular system. 
This task has been the primary focus of the Simulation Extensions Project, 
and, accordingly, has received the bulk of the available support. The 
Simulation Applications Task has studied issues related to generic surveil 
lance system performance drivers using tools developed under the Simulatio 
Enhancements Task. 

Section 2 summarizes the technical accomplishments of the Simulation 
Extensions Project and provides a guide to detailed techical documentation 
generated during the course of the project. 

I.2 Generic Requirements for the ROFT Program 

The development of space-based infrared surveillance systems and 
space-based lasers is subordinated primarily to the availability of very 
large, high precision optics. Current technology and existing facilities 
cannot produce in a timely fashion, high performance mirrors, both in qual 
ity and quantity, to respond to the needs of the systems envisioned. The 
Rapid Optics Fabrication Technology (ROFT) Program will attempt to discove 
and develop technologies which will significantly decrease the time it 
takes to produce optics. 

1 





The objective of this study reported in Section 3, is to provide 
broad requirements and constraints for the mechanical and thermal proper¬ 
ties of materials envisioned for the design and rapid fabrication of large 
space mirrors. By inspecting the constitutive equations relating mirror 
deformations to input disturbances, the relative importance of material 
parameters is identified. It is also recognized that the requirements on 
mirror material properties can vary widely depending on the space applica¬ 
tions. Surveillance systems, operating at cryogenic temperatures impose 
different demands on materials than do the laser systems for which the 
operational temperature is expected to vary between 300° K and 400° K. 


Dynamic environments exhibit significant differences as well, and this 
difference may impact the choice of other properties required of the mirror 
materials. 

It must be emphasized that this study does not presuppose a specific 
mirror concept with defined structural and configurational characteris¬ 
tics. While this approach is intentional for the findings to remain gen¬ 
eral, it limits our ability to obtain specific requirements for the contem¬ 
plated mirror materials. Furthermore the effectiveness of possible con¬ 
trols (whether thermal or mechanical), in terms of relieving demands on 
material properties can only be evaluated in the context of a defined mir¬ 
ror concept. As a result only estimates on mirror material requirements 
can be provided and they have been derived from simplified analytical ex¬ 
pressions and engineering experience in the design, modelling and analyses 
of space mirrors. 

For a large number of existing and potential materials, that include 
metals, composites and ceramics, all relevant properties have been tabu¬ 
lated for two extreme operational temperatures (150° K and 350“ K). The 
plots relating important material parameters can be useful in delineating 
domains of acceptability for a potential new material. The proposed mate¬ 
rial requirements can be viewed as somewhat flexible so as to permit trade¬ 
offs between properties for a particular design concept. Past experience 
has shown that the most likely candidate for such a mirror will probably be 
a thin face-sheet on actuators supported by a thermally stable substrate. 











However, it is hoped that dramatic changes and/or improvements in some 
thermal properties, especially, might lead to completely novel concepts 
that ultimately could be produced at the desired rapid rate. 


It should be remembered that the objective of the ROFT program is 
the rapid fabrication of high precision large space optics. The results of 
this study can only support this attempt by providing guidelines to the 
designer and fabrication engineer who selects possible materials. The gen¬ 
eric specifications of mirror materials, therefore, can be viewed as con¬ 
straints to the selection of a material amenable to rapid fabrication pro¬ 
cesses. Many other issues related to optical and perhaps electrical mate¬ 
rial properties need to be considered. Finally the fabrication process 
could be adversely affected by some mechanical and thermal properties of a 
material that is otherwise acceptable in the final product. 

1.3 Active Control of Space Structures 

Surveillance missions envisioned for future spaceborne optical 
systems strongly suggest the need for active control and embody several 
factors which—taken together—pose a formidable challenge from the 
perspective of control science and technology. First, the supporting 
structures for the optical system can be expected to be substantially non- 
rigid, exhibiting low and clustered characteristic frequencies of vibration 
and very small inherent damping. Moreover because of their anticipated 
geometric complexity, substantial deterministic uncertainty in available 
finite-element structural models will persist. Second, the disturbance 
environment may be severe, involving disturbances with nonnegligible power 
spectral density over a wide frequency range. Finally, high precision 
performance specification in regard to rapid reorientation and attitude 
regulation of the optical axis are certain. Such requirements imply the 
capacity for effective suppression of structural vibrations. 

The factors cited motivate the three main lines of effort in active 
control technology development that have been pursued under the present 
contract: (1) Synthesis of reduced-order control for vibration 










suppression; (2) System identification, with special attention to the 
problem of resolving closely-spaced modes; and (3) Generation of smooth 
large-angle-slew control. 

An overview of the work accomplished and principal contributions 
made in these three areas over the span of the contract is given in 
Section 4. The remaining sections detail the accomplishments of investi¬ 
gations carried out during the final semiannual reporting period. One of 
the principal efforts in reduced-order vibration control has been an 
experimental (not in hardware) approach to the study of mutual inter¬ 
actions between reduced-order modeling, active transducer (actuator and 
sensor) selection, and controller feedback strategy determination in the 
context of wideband disturbance accommodation for representative precision 
flexible structures. The specific insights obtained from this approach are 
discussed in Section 5; current results relate to suggested approaches for 
independent sensor selection, and to constraints imposed upon the synthesis 
process by the use of certain system identification techniques for 
reduced-order modeling. The enhancement of stability robustness to un¬ 
modeled high-frequency dynamics has been a major objective of the experi¬ 
mental synthesis study. Through the use of linear-quadratic-Gaussian (LQG) 
control employing loop transfer recovery (LTR), reduced order controllers 
have been obtained which stabilize the full 150-mode model of ACOSS Model 
No. 2. This successful application, and the limitations of the LQG/LTR 
approach in the present context, are presented in Section 6. An overview 
of our approach to the modeling of actuator dynamics, together with its 
connection to other relevant work, is given in Section 7. A principal 
focus of work in system identification has been the development of techni¬ 
ques for estimating the spectral characteristics of lumped-parameter 
systems with closely-spaced modes of vibration. Advances in this develop¬ 
ment are detailed in Section 8; special emphasis is given to the role of 
the sampling period used in obtaining input data for the identification 
process. Interest in identification processes for systems described by 
distributed-parameter models is rising because of current hardware develop¬ 
ments with distributed actuation and sensing devices for flexible structure 












vibration control. Theoretical development of 


efficient computational 


algorithm based upon cubic splines for parameter estimation in distributed 
systems is given in Section 9. Finally, some current advances in large 
angle slew control are discussed in Section 10. The particular topic 
emphasized is linear tracking controllers incorporating the capacity for 
accommodating discrete disturbances. An underlying thrust in this and 
prior work in large angle slew control has been to develop closed-form 
analytic solutions to the Riccati-like matrix equations arising in the 
mathematical development. This leads to considerable increases in the 
computational efficiency of related algorithms for computing the control. 
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SECTION 2 


SIMULATION EXTENSIONS 


Simulation Enhancements 

Over the period of performance of the ACOSS El 
ion Extensions Project has been focussed princ 
ulation Enhancements Task. The broad objectiv 
generalize those portions of the Draper Integ 


(DIS) which initially were too limited in scope or flexibility. Figure 
presents a simplified block diagram of the DIS. Within this category, t 
major effort has been directed toward freeing the DIS user from complete 
dependence on outside vendors for obtaining infrared (IR) background rad 
ance maps (i.e., "scenes") required for use in simulation studies. This 
effort, the Generic Scene Simulation (GENESSIS) subtask, has resulted in 
the augmentation of the DIS with an additional set of software modules a 
data bases which give the DIS user a limited but powerful capability for 
the generation and manipulation of synthetic terrestrial scenes. The 
GENESSIS models and software were developed by Photon Research Associate 
(PRA) under subcontract to CSDL. 

In addition to the GENESSIS work, the Simulation Enhancements Tas 
has supported enhancements to the Integrated Large Space Structures Simu 
t ion (ILS 3 ) as well as enhancements to the Signal Processor Simulation. 
As initially developed, the ILS 3 was mainly geared to address platform 
and structural issues arising in the context of a HALO-like system. The 
goal of the ILS 3 enhancements subtask was to develop a package of algo¬ 
rithms and software to enable the ILS 3 to analyze both multi-mirror ray- 
traced sensor systems as well as single-surface, paraboloidal or spheric 
reflector systems. 


The Signal Processor Simulation i 


developed originally 


of the DARPA-sponsored program at CSDL for end-to-end performance analys 
and simulation of the Mini-HALO system. During the course of this work 













Block diagram. 


















number of natural as well as potentially powerful generalizations of key 
algorithms suggested themselves. Some of these generalized algorithms have 
been implemented and studied as part of the Signal Processor Simulation 
Enhancements subtask. 

2.1.2 Simulation Applications 

Above and beyond the cases that have been examined as part of speci¬ 
fying, implementing and evaluating the enhancements to the simulation dis¬ 
cussed above, the simulation has been used to study several areas of gener¬ 
ic concern to space-based IR staring mosaic sensor systems. The most sig¬ 
nificant of these were the following: clutter leakage as a function of 
spectral interval and scene local time-of-day through the signal processor 
of a generic space-based IR staring mosaic sensor; and effectiveness as 
well as limitations of a class of algorithms proposed by D. Fried for the 
post-detection reduction of line-of-sight (LOS) jitter effects. 

2.2 Simulation Enhancements 

2.2.1 Generic Scene Simulation (GENESSIS) 

In support of the Draper Integrated Simulations Extensions effort, 
Photon Research Associates (PRA) has developed (under subcontract to CSDL) 
a computer code capable of generating and manipulating terrestrial scenes 
as a function of major surveillance system and mission parameters. This 
code (called GENESSIS) has the capability to interface with the Defense 
Mapping Agency (DMA) data base of terrestrial scenes as the source of scene 
input data. Consequently, the code is able to simulate any scene for which 
DMA data exists. 

The GENESSIS scene simulation is based upon a point-by-point algo¬ 
rithm, a single cycle of which consists of collecting (and in some in¬ 
stances, computing) inputs specific to a single point on the scene, calcu¬ 
lating the apparent radiance of that point from the collected inputs, and 
finally weighting and assigning the calculated radiance to the appropriate 
pixel in the observer's field-of-view. If the density of points is large 
enough, the scene will be properly sampled and the radiances computed by 
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repeated point calculations can be combined to produce an accurate pixel 
radiance map of the scene. The parameters of these radiance grid points 
are computed from the three-dimensional scene itself. 

Scene data consists of discrete altitude, material type pairs speci¬ 
fied at regular intervals on a planar rectangular grid. Continuous sur¬ 
faces are produced from the discrete scene data using a bi-cubic spline 
fitting technique. Point data can be computed from these surfaces at any 
desired spatial resolution. 

The computed apparent radiances consist of four terms combined add- 
itively. These are reflected solar, thermal emission, reflected skyshine 
and path radiance. Each major calculational operation is performed with a 
separate software module. The atmospheric, geometric and radiance modules 
have stand-alone capabilities, but are normally executed in sequence to 
produce a final result. 

The simulations' primary output is an N * M viewer-perspective pixel 
apparent radiance map. Diagnostic output is also available to check proper 
code execution. Figure 2-2 presents a block diagram of the GENESSIS soft¬ 
ware architecture. 

The GENESSIS code is comprised of six (6) modules (subroutine pack¬ 
ages) each with a single specific task. These are geometric, atmospheric, 
heat transfer, radiance, image and ephemeris. These are combined into 
three major modules, each with stand-alone capabilities. Modular stand¬ 
alone capabilities allow flexibility of operation while maintaining a sim¬ 
plicity of structure, user interaction and memory requirements. 

The geometric module performs shadowing and the viewer perspective 
projection of the scene. Its output is required by the radiance module. 

The atmospheric module supplies atmospheric parameters required by 
the radiance module. It is run least often since its output covers a wide 
range of solar and observer geometries. 
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The radiance module produces a viewer perspective, pixel apparent 
adiance map from information supplied by the atmospheric and geome- 
dules. It calls upon the heat transfer and image modules to pro- 
espectively, surface temperatures and viewer image. The heat trans 
ule currently does not have a stand-alone capability. Both the geo 
and radiance modules utilize the ephemeris module. 

User inputs are categorized according to purpose. These are geome- 
ensor and atmospheric. The elements of these are: 

A. Geometric Inputs 

1. The date and time of the simulation, used to compute the 
position of the sun, 

2. The latitude and longitude of the viewer subsatellite point 

3. The observer altitude in kilometers. 

B. Sensor Inputs 

1. The vertical and horizontal angular field-of-view, 

2. Focal plane rotation in degrees, 

3. The vertical and horizontal spatial resolution in meters. 

C. Atmospheric Inputs 

1. Atmospheric model (six LOWTRAN standard atmospheres), 

2. Aerosol model, 

3. Haze model, 

4. Visibility in kilometers. 

Scene data inputs consist of the altitude, material type pairs plus 
, atmospheric and reflectance data. The elements of these inputs 


1. Material thermal properties (solar absorptance, thermal emit- 
tance, thermal conductance and thermal mass) required by the 
heat transfer module, 
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2. Material in-band diffuse reflectance, 

3. Cloud in-band bi-directional reflectance, 

4. Surface level atmospheric properties (temperature, wind speed 
and humidity). 

Table 2-1 presents a summary of GENESS1S features. More detailed 
information on GENESSIS may be found elsewhere.1^,5] 

2.2.2 Enhancements to the Integrated Large Space Structures Simulation 

(1LS 3 ) 

The ILS^ structures/optics simulation is a tool for use in sim¬ 
ulating the performance and response of a flexible optical system subject 
to onboard and environmental disturbances and control system forces. The 
principal output of this program is a time-history of the system optical 
point-spread function (PSF) which can then be used as an input to the 
SBS^ sensor simulation. The efforts summarized here represent the first 
in a series of enhancements whose objectives were to broaden the types of 
systems which could be simulated by removing the modelling restrictions 
that were part of the original program. The overall structure of ILS^ is 
shown in Figure 2-3. It consists of three groups of basic functional mod¬ 
ules. In the first group, the node-point displacements of a finite-element 
model subjected to static or dynamic loads are computed. In the second 
group, these node-point displacements are transformed into optical wave- 
front errors. Finally, in the third group the wavefront errors are pro¬ 
cessed to compute LOS errors, PSF, etc. 


The focus of this effort has been the enhancement of the second 
group of modules, known as the interface programs. In these programs the 
previously computed displacements at the node-points of the finite-element 
model of the surface are transformed into optical wavefront errors at a 
rectangular grid of points. This is accomplished by interpolating the 
known node-point displacements to find the displacements at the points in 
the output grid for each reflective surface in the system. These displace¬ 
ments can then be combined using the information obtained from a raytrace 
of the system to determine the wavefront error at each point in the grid. 

1 3 








MAJOR USER-SPECIFIED VARIABLES 

- Observer altitude, zenith angle 

- Spectral interval (2.5-13.0 Mm) 

- Atmospheric model (LOWTRAN family) 

- Field-of-view location (within overall scene) 

- Scene angular resolution 

- Continuously variable surface level wind speeds and 

temperatures 


FIVE REPRESENTATIVE TERRAIN DATA BASES 

- California Coast (north of Santa Cruz) 

- Brooks Range Mountains of Alaska 

- Arctic Tundra 

- Middle East 

- Central Europe 


FOUR CLOUD SCENES 

- Stratus 

- Altocumulus 

- Cirrus 

- Cumulus 


EACH DATA BASE TOPOLOGY EXTRACTED FROM DEFENSE MAPPING AGENCY 
(DMA) DATA BASE 


EACH DATA BASE REGISTERED WITH LANDSAT DATA TO OBTAIN MATERIALS 
MATERIALS ASSIGNMENTS 

- Limited to 14 material types 


TERRAIN SCENE SIZE: APPROXIMATELY 40 km * 40 km 



























The principal result of this effort has been to generalize the ILS 3 so that 
if the system contains only one mirror (which is spherical or paraboloidal) 
and a point feed, the wavefront errors for on-axis operation can be comput¬ 
ed directly from the geometry of the system without the need for a ray- 
trace . 


The methods used to perform the interpolation of the node-point dis¬ 
placements used to find the displacements on the optical grid are described 
in detail elsewhere.The interpolation actually involves two steps: 

1) a search to determine where on the finite-element model each optical 
grid point lies; and 2) the computation of the interpolation coefficients 
for each point. The algorithms used in both steps have been designed to 
allow any type of finite-element model for the mirror surface. 

2.2.3 Enhancements to the Signal Processor Simulation 

The signal processor simulation is a key part of the Draper Inte¬ 
grated Simulations (DIS). The DIS is a sophisticated analysis tool dev¬ 
eloped under DARPA support for overall evaluation and performance assess¬ 
ment of space-based surveillance systems. The DIS models in detail the 
mechanical, optical, control, signal collection and signal processing sub¬ 
systems in a highly interactive fashion. It is composed of a number of 
large data bases (including a library of simulated earth scenes) and two 
principal simulation packages: the Integrated Large Space Structures Sim¬ 
ulation (ILS 3 ); and the Space-Based Surveillance Systems Simulation (SBS 3 ). 
ILS 3 simulate the space platform; it assesses the geometrical and pointing 
performance of satellites and their sensors as influenced by thermal, 
structural, dynamic, and control forces. SBS 3 uses the ILS 3 outputs and 
models the performance of space-based surveillance systems with respect to 
their intended missions using nominal and perturbed values for the target, 
background, optics, LOS, optical filters, focal plane, and signal process¬ 
or. Figure 2-4 presents a simplified block diagram of the DIS. 

Initial development of the Signal Processor Simulation was driven by 
a requirement from an earlier DARPA-sponsored program at CSDL to model the 
HALO Signal Processor (HSP). As indicated in Figure 2-5, the HSP was com- 
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High Resolution Acquisition 
Processor (HRAP) 


Acquire Point-Source Targets 


High Resolution Track Processor 
(HRTP) 


Track Point-Source Targets 
Acquired by HRAP 


Low Resolution Acquisition 
and Track Processor (LRTP) 


Acquire and Track Resolved 
Targets (e.g. Plumes) 


Tactical Processor (TACPRO) 


Detect Explosions and Other 
Tactical Events 


Figure 2-5. Components of Draper Signal Processor Simulation (DSPS) 
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prised of four separate processor modules. Subsequent work under the Sim¬ 
ulation Extensions Project has focussed on the High Resolution Acquisition 
Processor (HRAP) and the Low Resolution Acquisition and Track Processor 
(LRTP), and for efficiency, these two simulations have been combined into 
one software package. The simulation user has the option of selecting the 
HRAP or the LRTP for processing a given set of input data frames (the HRAP 
and LRTP run independently, not interactively). In addition, there are 
numerous other options available to the user, comprising a menu which 
specifies the signal processor configuration for a given run. In fact, by 
tailoring the input parameters, the user can implement a number of differ¬ 
ent signal processor configurations which go far beyond either the HRAP or 
the LRTP. The CSDL HRAP/LRTP simulation is divided into a number of mod¬ 
ules, which generally perform operations on the input data in a sequential 
fashion, frame by frame. Depending upon the processing sequence specified 
by the user, some of the modules may or may not be used. 
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The first of the modules performs either temporal integration of 
successive frames or spatial integration within each frame (or both). 


The next module is the differencing module, which computes the 
zeroth-, first-, second-, or third-order differences (as specified by the 
simulation user) for a set of input data frames. A difference is essen¬ 
tially a derivative of the signal amplitude with time; zeroth order cor¬ 
responds to raw data, first-order differencing corresponds to a first deri¬ 
vative, and so on. 

Next is the thresholding module, which detects potential targets by 
comparing one or more difference frames with templates (normalized to a 
predetermined threshold level) that correspond to the minimum expected 
pattern in the differencing that a moving target would generate. Varia¬ 
tions on this technique include simple thresholding (which employs only one 
frame at a time), Multi-Threshold Detection (or MTD, which employs a multi- 
-frame template^, and Rationalized Exceedance Thresholding (RET, which is 
similar to MTD but has more rigorous requirements). A variety of thresh¬ 
olding templates can be used, depending on the difference order employed. 
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The threshold level setting module is called by the thresholding 
module in order to determine the signal level to which the thresholding 
templates will be normalized. The threshold level can be read in, computed 
on the basis of a noise estimate, or set to be proportional to the standard 
deviation of the current frame of data or differencing. The threshold 
level can be periodically reset by a Constant False Alarm Rate (CFAR) loop. 

Clustering and centroiding follow threshold level setting and thres¬ 
holding. The clustering module selects groups of thresholded hits for 
centroiding. The centroiding module either computes the moment of each 
group of hits or will select just the peak value. There are several 
options for the amplitude basis for clustering and centroiding, or, if the 
user chooses, these modules can be bypassed altogether. 

The output of the centroiding module is fed to the acquisition mod¬ 
ule. The acquisition module can employ the two tier algorithm or the one 
tier algorithm for high resolution or low resolution targets, respectively. 
The acquisition algorithms form tracks by associating any set of sequential 
and contiguous hits subject to certain rules. There are flexible criteria 
for determining when a track has a sufficient number of data points to be 
considered acquired. 

Tracks which are generated by the acquisition module are handed over 
to the least-squares line fitting module or the correlation coefficient 
module. These modules eliminate false tracks by rejecting any tracks which 
have, 1) points at a greater than acceptable distance from the track's 
least-squares generated line, or 2) a correlation coefficient less than 
some minimum allowable limit. 

Table 2-2 summarizes the major user variables for the enhanced sig¬ 
nal processor simulation. Figure 2-6 illustrates the collection of modules 
which comprise the CSDL HRAP/LRTP simulation. Further information can be 




-2. Major user variables for the enhanced signal proce 
simulation 


VARIABLE 

RANGE 

Difference Order 

• Zeroth (i.e., no frame differencing) 

• First order 

• Second order 

• Third order 

Threshold Level Selection 

• User specified 

• (Noise estimate) * (Threshold-to- 
noise Ratio) 

Track Acquisition Algorithm 

• Low resolution type 

• High resolution type 

Difference order for use in 
low resolution noise estimate 

• First order 

• Second order 

• Third order 

Thresholding Polarity 

• Negative contrast only 

• Positive contrast only 

• Arbitrary contrast 

Thresholding Algorithm 

• Simple 

• HSP Multi-Threshold Detection (MTD) 

• CSDL MTD 

• Idealized MTD 

• Rationalized exceedance thresholding 

Constant False Alarm Rate 
(CFAR) Option 

• Number of allowed false alarms and 
maximum number of iterations 

Centroid Algorithm 

• Conventional 

• Peak value 

Cluster/Centroid Basis 

• Raw data 

• Difference data 

• Threshold exceedance 

System Track Criterion 

• Ml hits out of N1 frames (Ml < Nl) 

Acquired Track Criterion 

• M2 hits out of N2 frames (M2 < N2) 

Frame Temporal Integration 

• Number of frames to be integrated 

Frame Spatial Integration 

• Number of pixels to be integrated 

Noise Estimation 

• Number of frames to be used 

False Track Rejection 

Algorithm 

• Least-Squares 

• Correlation coefficient 

Acquired Track Linearity 

• Maximum allowed deviation from 
best-fit straight line 
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Figure 2-6. HRAP/LRTP Modules and data flow. 


















2.3 Simulation Applications 
2.3.1 Clutter Leakage Study 

The Simulation Enhancements described above have considerably ex¬ 
panded the scope and capability of the DIS to conduct end-to-end perform¬ 
ance assessments of space-based staring infrared surveillance systems. As 
an illustration of this new capability, the DIS has been used to study the 
effect of choice of spectral interval and selection of scene local time of 
day on clutter leakage through a generic signal processor.The basic 
approach was to begin with the California coast (Santa Cruz) data base and 
to use GENESSIS to produce two scenes in the 3.6-4.0 pm band and two scenes 
in the 8.0-9.0 Pm band, each pair of scenes corresponding to 8 AM and to 
noon local time, respectively. A nominal surveillance system mission and 
configuration were assumed, and the above scenes were processed through the 
DIS under varying levels of (LOS) jitter. Tables 2-3 and 2-4 summarize the 
scenario, sensor and signal processor configurations for this investiga¬ 
tion. These parameter values reflect the basic configuration for a system 
whose principal objective is the detection of strategic aircraft against a 
terrestrial background. The figure of merit used to measure clutter leak¬ 
age in the signal processor was average number of threshold exceedances per 
frame. 


Since the average number of threshold exceedances per frame is a 
function of threshold level in the signal processor as well as the level of 
LOS jitter, the results of this study may be presented from those two 
points of view. Figure 2-7 illustrates the dependence on LOS jitter level 
of the average number of threshold exceedances per frame. The particular 
threshold level selected is representative of what would be used to acquire 
a target with a relatively high signal-to-noise ratio. The general trend 
of the data is such that, for a given level of LOS jitter, higher numbers 
of exceedances are found in the LWIR band (i.e., 8.0-9.0 pm) than in the 
SWIR band (i.e., 3.6-4.0 Pm). In either band, the higher numbers of 
exceedances always occur at noon local time as opposed to 8 AM. 
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Table 2-3. Scenario and sensor specificat 
Santa Cruz investigation. 


Target 

None 

Backgrounds 


-Location 

Santa Cruz (GENESSIS) 

-Spectral Bands 

3.6-4.0 pm 

8.0-9.0 pm 

-Times of Day (Local) 

8 AM 

NOON 

Spacecraft 


-Drift 

None 

-Jitter 

0.1 prad rms 

0.5 Prad rms 

1.0 prad rms 

-Altitude 

35,700 km 

-Boresight 

NADIR 

Optics 


-Aperture Diameter 

1.4 m 

-Obscuration Ratio 

0.35 

-Thruput 

0.135 

-Point-spread Function 

Gaussian (o = 5.01 Prad) 

Focal Plane 


-Detector Geometry 

25 x 25 Array of 100 pm 
(20 prad) Square Pixels 

-Frame Duration 

3.4 sec 

-D* 

2.8 (10) 12 cm /Hz/W 

-Responsivity non-uniformity 


.Fixed Pattern 

OX 

.Random 

15% 
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Difference Filter 

3rd Order 

Threshold Level 

User Input 

Assumed Target Contrast 

Negative 

Threshold Templates 

CSDL MTD 

Threshold Algorithm 

CSDL Rationalized 
Exceedance Thres¬ 
holding (Version 1) 

Cluster Basis 

Threshold Exceedances 

Centroid Basis 

Threshold Exceedances 

System Track Criterion 

3 out of 3 

Acquired Track Criterion 

10 out of 12 

System Track 

Threshold Scaling 

0.5 

Least Squares Track 
Linearity Criterion 

±1.8 Pixels 



















Figure 2-8 shows the typical dependence of average number of thresh¬ 
old exceedances per frame on threshold level for a given value of LOS jit¬ 
ter. Many of the same trends discussed above are also evident in this 
data. Decreasing the threshold results in a rapid increase in threshold 
exceedances. At a given threshold value, higher numbers of threshold ex¬ 
ceedances always occur in the LWIR band as opposed to the SHIR, and, within 
each band, higher numbers of exceedances occur at noon in comparison with 8 
AM. Increasing the LOS jitter level results in a translation upward of the 
curves at higher threshold values. In addition, there is a compression 
within each spectral band of the exceedance curves corresponding to noon 
and to 8 AM. 

2.3.2 Evaluation of Algorithms for the Post-Detection Compensation Of 

Line-of-sight (LOS) Jitter Effects 

The Draper Integrated Simulations (DIS) in general provides a nat¬ 
ural test bed for examining the behavior of current and proposed space- 
based surveillance systems. In particular, since LOS stability has been 
found to be a major performance-limiting factor for a broad class of 
downward-staring surveillance systems, the DIS has been used to evaluate a 
technique proposed by D. Fried^l f or reduction of the LOS jitter 
effects. 

Moving targets are most commonly detected by mosaic sensors through 
differencing of successive frames. Unfortunately, the changes in the de¬ 
tectors' output are not only due to moving or changing targets, but also to 
jitter and drift of the sensor LOS as well as to detector noise. Under 
certain assumptions, there is some retrievable information about the jitter 
contained in the sequence of frame outputs from the detector array. Thus, 
assume that the focal plane data rate is higher than the cut-off frequency 
of the jitter and that the jitter amplitude is much smaller than the 
spatial resolution of the sensor optics/detector. Then, in the absence of 
detector noise, the movement of the sensor between two arbitrary frames can 
be estimated and this estimate used to correct the output before further 
processing to extract the target. 





The estimation method used is essentially that proposed by D. Fried, 
et. al. Figure 2-9 is a schematic of a portion of the mosaic array which 
has moved by 6x end <$y between two frames at times t and t' respectively. 

If S A ^respectively, S A t ,) and Sg ^respectively, Sg t ,) are the voltage 
outputs of adjacent detectors A and B at time t (respectively, t 1 ), and A 
is the center-to-center detector spacing, then, with the restrictions 
mentioned previously, 


Using these approximations, it can be shown that 


\ / 


is a good estimate of the displacements 6x. Tie symbols < > represent the 
averaging operator over a portion of, or perhaps the entire, focal plane. 

A similar expression can be obtained for Ay by using pairs of detectors A 
and C. With the same restrictions as before, it follows that 
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Depending on the scenario selected, there are several ways of imple¬ 
menting the correction of the jitter effects. One method is to employ a 
stabilization mirror with a feedback using the estimated values of the 
jitter, but this approach may be difficult to implement in the case of sen¬ 
sors with a wide field-of-view. The method selected here accomplishes the 
correction through data processing. It has the advantage of allowing for 
localized correction of the jitter effects when the field-of-view is large 
and/or the aperture is large and requires segmentation of the optical 
elements. 

In order to make the correction scheme possible, some of the hard¬ 
ware parameters of the focal plane must be appropriately selected. For the 
estimation of the jitter, the information regarding the spatial slope of 
the pixel ouput as a function of position must be preserved. Therefore, if 
the focal plane has a reset function, it must be suppressed. This could be 
achieved by selecting an initial reference voltage approximately equal to 
the mid-value over the array of the voltage outputs of the multiplexer for 
the scene considered and by selecting an off-chip amplifier gain small 
enough so that no detectors reach saturation. Note that this second point 
creates the need for a finer discretization level in the A/D conversion 
process. 

A study of the influence of the PSF size has shown that a larger PSF 
allows for better improvement in clutter reduction compared to a small PSF. 
Also, for highly structured backgrounds it was found that the fixed pattern 
noise induces a reduction of the estimated value of the jitter, and there¬ 
fore, a poorer correction when compared to a homogeneous background. For 
the cases considered, the jitter correction scheme was used over a range of 
rms jitter values from 0.1 to 1 microradians. In this range it was found 
that, for highly structured backgrounds and an optimized PSF size, the 
effect of the correction scheme was equivalent to a reduction of the LOS 
jitter by approximately a factor of five. 
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SECTION 3 

GENERIC MATERIAL REQUIREMENTS FOR THE ROFT PROGRAM 

3.I Introduction 

The development of space-based infrared surveillance systems and 
space-based lasers is subordinated primarily to the availability of very 
large, high precision optics. Current technology and existing facilities 
cannot produce in a timely fashion high performance mirrors, both in qual¬ 
ity and quantity, to respond to the needs of the systems envisioned. The 
Rapid Optics Fabrication Technology (ROFT) Program will attempt to discover 
and develop technologies which will significantly decrease the time it 
takes to make such optics. The goal of the ROFT program is to make possi¬ 
ble a significant increase in mirror production rate, from the traditional 
one square meter per year to about 100 M^ per year. Such an increase 
represents a giant step forward and it will be necessary to exhaust all 
possible avenues capable of enhancing rapid fabrication procedures of high 
quality optics. 

The technology has been developed to make single mirrors having 
acceptable weight and appropriate performance characteristics for laser and 
multi-mission surveillance systems. Programs like HALO and LAMP recently 
have demonstrated the ability to make mirrors capable of controlling their 
figure to the required surface accuracy under the disturbances of both 
orbital and control environments. The design of these mirrors was pre¬ 
occupied mainly with adequate performance in orbit, but did not address the 
question of making these high quality optics in large quantities. Signif¬ 
icant departure in the design approach and in the fabrication procedures of 
these mirrors must be incorporated. In attempting to achieve the desired 
high rate of production, new technology must be concerned with each step of 
the mirror fabrication process: material production, blank fabrication, 
mirror generation, grinding, polishing, metrology and assembly. In order 
to remain completely general in its approach, the ROFT program will also 
address unconventional system concepts which could have the potential to 
circumvent the need for large high quality traditional mirrors. Addition- 
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ally, some effort will be channeled towards the rapid development of 
optical metering trusses since the fabrication here is also crucial, but 
not nearly as much as for the mirrors. 

The overall objective of this CSDL effort is to provide broad 

requirements and constraints for the mechanical, thermal and optical prof 

erties of materials envisioned for the design and fabrication of new opti 
concepts. This study, which is a subset of the ROFT program plan, provic 
only generic specifications so as not to restrict the freedom of the de¬ 
signer in his attempt to develop a new technology and thus enhance the 
chances for the successful development of a rapid fabrication process. 

In addressing the design of the space mirrors, two distinct applic 
tions have been considered: (a) optics appropriate for use as the main 
optics in an infrared (IR) surveillance systems; (b) optics appropriate I 
use as the primary mirror of a beam expander for a space-based laser sys¬ 
tem. The incident flux in the aperture plane of the IR system differs 

greatly from that of a laser and the respective operating temperatures 
(150° K for IR and 300° to 400° for lasers) pose different problems for 
each system. 

The general requirements for mirror materials for both applicatioi 
are that the mirror surface be of good optical properties. These materii 
should have a very low coefficient of thermal expansion and a very small 
variation of this coefficient. They should exhibit a low density combine 
with a high modulus of elasticity and micro-yield strength. A very impoi 
tant requirement for the ROFT program is that a proposed mirror material 
must be amenable to rapid fabrication or replication at low cost and rem< 
within reasonable weight constraints. Finally the mirror concepts must I 
deployable in space and exhibit very little outgassing. 

This section summarizes the approach used in this study to define 
the mirror material requirements and presents some of the important re¬ 
sults. Most of the technical derivations have been presented in some 
detail earlier and the interested reader should refer to Section 3 of th< 
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ACOSS Eleven Fifth Semiannual Technical Report.* The conclusions and 
recommendations presented here include additional design suggestions that 
could be explored with the promise of significant pay-offs. 


3.2 Materials Characterizatii 


Approach and Results Summary 


The performance of a large high precision space mirror is subordin¬ 
ated to a number of variables that are, in some cases, closely interre¬ 
lated. Starting with the definition of a structural concept, fabrication 
and assembly procedures must be devised so as to permit the rapid produc¬ 
tion of the mirror. In addition the material selected must exhibit opti¬ 
cal, elastic and thermal properties that will insure an acceptable perform¬ 
ance for the mirror exposed to the rigors of environmental and control dis- 


The performance of a mirror structure, regardless of its structural 
concept, can be estimated from the equations expressing the distortions as 
a function of the structural characteristics and the expected loadings. 

The mirror deformation d is affected by the three following contributions: 


(1) geometry and boundary conditions of the structure (Fg), 

(2) the material mechanical and thermal properties (F^), and 

(3) the loading conditions (Fl) in an expression of the form 


The approach suggested here attempts to uncover those materials 
properties (entering the equation through F^) that affect most adversely 
a mirror exposed to the various loading conditions expected in space appli¬ 
cations. 

The behavior of a mirror was, for simplicity, associated with the 
central deflection or sag d of a thin circular plate with an initial curva- 

* ACOSS Eleven Fifth Semiannual Technical Report, CSDL-R-1687, March 1984. 










Cure of radius R. The plate is assumed simply supported along its periph¬ 
ery. (See Figure 3-1) The equations can be expressed either in terms of 
constant mass per unit area (m) or constant mirror thickness (h). 



Figure 3-1. Mirror description. 

The equations expressing the sag of a typical mirror plate have been 
derived for the following loading conditions: 

• uniform pressure (examples are gravity and polishing loads) 

• harmonic forcing function (dynamic responses) 

• thermal steady state and transient responses 

• thermal conditions with a inhomogeneities and laser beam spikes 

• bulk temperature changes due to orbit insertion. 

The derived equations furnish the engineer with a useful tool to 
assess the respective merits of various mirror materials. It is apparent 
that materials properties such as a, Aa, k, C, E and P determine to a 
great extent the mirror distortions likely to be encountered in space 
applications. A large number of candidate materials were examined that 
include: (a) metals, (b) composites, and (c) glassy materials. Some of 









these materials are well known while others are still developmental. Table 
3-1 presents the list of 34 materials with their elastic and thermal prop¬ 
erties at the 150° K expected in IR applications. 


Table 3-2 was prepared as an attempt to place some desirable bounds 
on the most important material properties and critical parameters involving 
several constants. They are best estimates based on both analytical de¬ 
rivations and engineering experience with space mirror designs and analy¬ 
ses. It is understood that these constraints are not mandatory but can 
serve as reasonable guidelines to initiate trade-offs between them. The 
table also identifies the existing materials that exhibit the best proper¬ 
ties in each category 


In order to appreciate where the existing material properties fall 
with respect to desired values, several plots relating critical parameters 
have been generated. It is possible to see more directly which materials 
proposed for mirror applications fall within a region of acceptability. 

For example, Figure 3-2 identifies the materials on the left of the 

p(l_ v 2) a _15 

-jjj- • = 10 line as the most promising for a IR surveillance sys¬ 
tem. Other similar plots presented earlier for other important parameters 
provide a more direct insight into the merits of candidate mirror materials 
and identify interesting and promising groups. 


3.3 


Conclusions and Recommendations 


The objective of this study was to provide broad requirements and 
constraints for the mechanical and thermal properties of materials envi¬ 
sioned for the design and rapid fabrication of large space mirrors. By 
inspecting the constitutive equations relating mirror deformations to input 
disturbances, the relative importance of material parameters was identi¬ 
fied. It was also recognized that the requirements on mirror material 
properties can vary widely depending on the space applications. Surveil¬ 
lance systems, operating at cryogenic temperatures impose different demands 
on materials than do the laser systems for which the operational tempera¬ 
ture is expected to vary from 300° K to 400° K. Dynamic environments ex- 
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Desirable elastic and thermal properties for potential mirror material 






















hibit significant differences as well, and this difference may impact the 
choice of other properties on the mirror materials. 

It must be emphasized that this study did not presuppose a specific 
mirror concept with defined structural and configurational characteris¬ 
tics. While this approach was intentional for the findings to remain gen¬ 
eral, it limits our ability to obtain specific requirements for the contem¬ 
plated mirror materials. Furthermore the effectiveness of possible con¬ 
trols (whether thermal or mechanical), in terms of relieving or modifying 
demands on material properties can only be evaluated in the context of a 
defined mirror concept. As a result only estimates on mirror material 
requirements could be provided and they have been derived from simplified 
analytical expressions and engineering experience in the design, modelling 
and analyses of space mirrors. 

For a large number of existing and potential materials, that include 
metals, composites, glass and ceramics, all relevant properties have been 
tabulated for two extreme operational temperatures (150* K and 350“ K). 

The plots relating important material parameters can be useful in delinea¬ 
ting domains of acceptability for a potential new material. One should 
emphasize that the proposed material requirements should be viewed as some¬ 
what flexible so as to permit tradeoffs between properties for a particular 
design concept. Past experience has shown that the most likely candidate 
for such a mirror will probably be a thin face-sheet on actuators supported 
by a thermally stable substrate. However, it is hoped that dramatic 
changes and/or improvements in some thermal properties, especially, might 
lead to completely novel concepts that ultimately could be produced at the 
desired rapid rate. 

A promising approach would be one that recognizes in the mirror 
design and its responses the intricate relationships between the elastic 
and thermal characteristics of the structural materials. The structural 
concept and configuration acquires importance as well since they could be 
manipulated in such a way as to optimize the thermo-elastic responses of 








the mirror. Some ideas, as the ones described below, could be explored to 
assess their potential payoffs. 


(1) Take advantage of the use of different materials where dii 


similar thermal and elastic properties combine to enhance 


structural behavior and optical performance. 


(2) Based on an optimal range for material properties (for a 

specific mirror concept) modify properties of constitutive 
materials of a nominal design to match desired ideal mix of 
parameters. (Today's materials technology shows promise in 
tempering with "nominal" material characteristics.) 


(3) Investigate the many possible thermal control capabilities that 
would advantageously modify mirror temperature distributions 
and ultimately enhance optical performance. 


It should be remembered that the objective of the ROFT program is 
the rapid fabrication of high precision large space optics. The results of 
this study can only support this attempt by providing guidelines to the 
designer and fabrication engineer who selects possible materials. The gen¬ 
eric specifications on mirror materials, therefore, can be viewed as con¬ 
straints to the selection of a material amenable to rapid fabrication pro¬ 
cesses. Many other issues related to optical and perhaps electrical mate¬ 
rial properties need to be considered. Finally the fabrication process 
could be adversely affected by some mechanical and thermal properties of a 
material that is otherwise acceptable in the final product. 
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SECTION 4 


ACTIVE CONTROL OF FLEXIBLE SPACE STRUCTURES: OVERVIEW 


4.1 Scope 

Surveillance missions envisioned for future spaceborne optical 
systems strongly suggest that active control will be needed to ensure 
that the mission objectives are achieved. Moreover, such missions embody 
several factors which—taken together—pose a formidable challenge from 
the perspective of control science and technology. First, the supporting 
structures for the optical system can be expected to be substantially 
nonrigid, exhibiting low and clustered characteristic frequencies of 
vibration and very small inherent damping. Moreover because of the anti¬ 
cipated geometric complexity of such structures, substantial determinis¬ 
tic uncertainty in finite-element models of the structures will persist. 
Second, the disturbance environment may be severe, involving disturbances 
with nonnegligible power spectral density over a wide frequency range. 

For example, the dynamics of fluid flow associated with the cooling of 
optical surfaces may produce such disturbances. Finally, performance 
specifications requiring high precision in regard to rapid reorientation 
of and attitude regulation of the optical axis are certain. The ability 
to satisfy these performance specifications requires an (implied) capa¬ 
city for effective suppression of structural vibrations. 

The factors cited above motivate the three main lines of effort in 
active control technology development that have been pursued under the 
present contract: (1) Synthesis of reduced-order control for vibration 











suppression; (2) System identification, with special attention to the 
problem of resolving closely-spaced modes; and (3) Generation of smooth 
large-angle-slew control. 

4.2 Major Accomplishments 

Significant results have been obtained in the course of pursuing 
each of the main lines of effort. They are summarized in turn. 

4.2.1 Vibration Control 

The principal objective in this effort has been to synthesize 
active controllers that maintain specified optical performance (including 
full-order stability) in the presence of wideband disturbances. Initial 
experience Indicated the need for a systematic examination of the mutual 
interactions between three principal elements of the overall synthesis 
process: reduced-order modeling; active transducer (actuator and sensor) 
selection; and controller feedback strategy determination [Ref. 4-3 
(Sec. 5)]. A new analytical approach to the selection of active trans¬ 
ducers from among a large class of candidates was conceived [Ref. 4-2 
(Sec. 3)]. The approach was then developed Into an efficient tool for 
examining the influence of active transducer placement (Including non- 
colocated configurations) upon closed-loop stability and performance 
[Refs. 4-4 (Secs. 2, 3), 4-5 (Sec. 5)]. An encouraging degree of 
stability robustness to unmodeled high-frequency dynamics, attributable 
to active transducer placement alone, was demonstrated. Finally, by 
applying recent insights |Ref. 4-6] to modifying the controller feedback 
strategy, closed-loop stability of the full-order (150-mode) model of a 
spaceborne optical system [Refs. 4-7, 4-8] with a reduced-order (18-mode) 
controller was demonstrated (Section 6 herein). During the latter demon¬ 
stration, the inadequacy of current technology for maintaining perfor¬ 
mance robustness was clarified. 

Several notable theoretical contributions relevant to vibration 
control have been made independent of the synthesis demonstration 








studies. New insights into the modeling of actuator dynamic processes 
[Ref. 4-3 (Sec. 3), and Section 7 herein] and sampled-data processes 
[Refs. 4-4 (Sec. 4), 4-5 (Sec. 6)] have been generated. The key to these 
insights is the use of natural second-order dynamic models—which retain 
physical insight—rather than first-order state-space models—which 
obscure such insight. 

Although most of the theoretical and applied investigations under 
the contract employed lumped-parameter (finite-element) structural 
models, some significant theoretical results for systems retaining a 
distributed-parameter representation were obtained. Such investigations 
appear timely in view of the increasing level of interest in developing 
distributed hardware devices for the control and estimation of vibrating 
systems [e.g., Ref. 4-9]. Convergence of an algorithm employing a 
sequence of (finite-dimensional) stochastic optimal output feedback com¬ 
pensators for the control of a distributed system was established [Refs. 
4-1 (Sec. 7), 4-10]. 

4.2.2 System Identification 


Motivation for this effort lies in the unavoidable (deterministic) 
uncertainty in the available models for complex structures, which leads 
to the potential for instability under feedback control. The principal 
objective of efforts in this direction has been to develop estimation 
processes appropriate for generating reduced-order models of flexible 
space structures for the purpose of controller synthesis. A new approach 
to spectrum estimation for lumped-parameter systems based on discrete¬ 
time signal processing ideas was conceived [Ref. 4-3 (Sec. 2)]. The 
approach is specifically designed to accommodate the special characteris¬ 
tics exhibited by flexible space structures (e.g., closely-spaced vibra¬ 
tion nodes). It is particularly distinguished by the use of the sampling 
period (for input data to the identification process) as a parameter to 
aid in high resolution. The approach has been developed into a 
systematic and effective tool for high resolution spectrum estimation 
[Ref. 4-5 (Sec. 7) and Sec. 8 herein; cf. also Sec. 5.2.3 herein]. 
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Several approaches to parameter estimation in distributed- 
parameter systems have also been developed [Ref. 4-5 (Sec. 8), and Sec. 9 
herein)]. Convergence of solutions to finite-dimensional identification 
problems formulated with cubic-spline approximations has been proved, and 
numerical results have been generated to demonstrate a variety of appli¬ 
cations. 

4.2.3 Large Angle Slew Control 

The principal objective of this effort has been to synthesize 
active control strategies that enable large-angle spacecraft reorienta¬ 
tion between specifically constrained end conditions with a minimal exci¬ 
tation of structural vibrations. Although the general approach being 
developed is applicable to nonlinear problems [Ref. 4-11], the emphasis 
of work in this direction has been placed upon problems with linear 
dynamics (e.g., tracking) so as to obtain sharper results. In particu¬ 
lar, this allows the use of techniques in (closed-loop) linear multivari¬ 
able control [e.g., Ref. 4-12] to be employed, in addition to techniques 
in (open-loop) optimal control [e.g.. Ref. 4-13] applicable to nonlinear 
problems. Simulation results representing motion of a rotating rigid hub 
with four flexible appendages and demonstrating both open-loop and 
closed-loop control have been given [Refs. 4-1 (Sec. 8), 4-2 (Sec. 7), 

4-4 (Sec. 5)]. A significant aspect of the results relating to linear 
closed-loop slew control is the obtaining of closed-form solutions for 
various matrix Riccati-like equations arising in the mathematical 
development. Such representations enable considerable improvement in the 
efficiency of related algorithms used for computing the control. 


4.3 Summary of Current Results 

The remaining sections of the present report detail the accom¬ 
plishments of investigations carried out during the final semiannual 
reporting period. One of the principal efforts in reduced-order vibra¬ 
tion control has been an experimental (not in hardware) approach to the 
study of mutual interactions between reduced-order modeling, active 









transducer (actuator and sensor) selection, and controller feedback 
strategy determination in the context of wideband disturbance accommoda¬ 
tion for representative precision flexible structures. The specific in¬ 
sights obtained from this approach are discussed in Section 5; current 
results relate to suggested approaches for independent sensor selection, 
and to constraints imposed upon the synthesis process by the use of high- 
resolution spectrum estimation techniques for reduced-order modeling. 

The enhancement of stability robustness to unraodeled high-frequency 
dynamics has been a major objective of the experimental synthesis study. 
Through the use of linear-quadratic-Gaussian (LOG) control employing loop 
transfer recovery (LTR) [Ref. 4-6], reduced-order controllers have been 
obtained which stabilize the full 150-mode model of ACOSS Model No. 2. 
This successful application, and the limitations of the LQG/LTR approach 
in the flexible structure control context, especially in regard to per¬ 
formance robustness, are presented in Section 6. A mature perspective on 
our approach to the modeling of actuator dynamics, together with its 
connection to other relevant work, is given in Section 7. As indicated 
in the preceding subsection, a principal focus of work in system identi¬ 
fication has been the development of techniques for estimating the spec¬ 
tral characteristics of lumped-parameter systems with closely-spaced 
modes of vibration. Advances in this development are detailed in Section 
8; a thorough account of the parametric use of the data sampling period 
to improve spectral resolution is given. Theoretical development of a 
more efficient computational algorithm based upon cubic splines for para¬ 
meter estimation in distributed systems is given in Section 9. Finally, 
some current advances in large angle slew control are discussed in 
Section 10. The particular topic emphasized is linear tracking con¬ 
trollers incorporating the capacity for accommodating discrete distur¬ 
bances. Analytic representations of closed-loop solutions for matrix 
Riccati-like equations are discussed in detail. 









The emphasis in the work reported here has been upon Independent 
development—In some-depth—of the three main lines of effort Identified 
in Section 4.1. Results obtained reflect substantial maturity in each of 
the developments. However, for system-level applications of the techno¬ 
logy represented, an integrated approach embracing each of three lines of 
effort is required. Except for the analytical results of Sec. 5.2.3 of 
the present report, no such interdisciplinary work has been undertaken. 
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SECTION 5 


ACTIVE CONTROL SYNTHESIS FOR WIDEBAND DISTURBANCE ACCOMMODATION 
IN PRECISION FLEXIBLE SPACE STRUCTURES 

5.1 Motivation 

5.1.1 The Synthesis Problem 

Surveillance missions envisioned for future spaceborne optical 
systems lead quite directly to some challenging problems In structural 
control. Because of size, material properties, and certain structural 
symmetries, the support structure housing the optical train can be 
expected to exhibit low and clustered characteristic frequencies—low 
enough to be excited by any attitude controller—and very small inherent 
damping. To satisfy mission requirements, the spacecraft must have the 
capacity for rapid reorientation, timely suppression of structural 
vibrations, and regulation of the optical axis along a specified 
direction to high precision. Finally, this performance is to be 
maintained in the presence of disturbances having nonnegllglble power 
spectral density over a wide band of frequencies—for example, 
disturbances generated by fluid flow for mirror thermal regulation. 
Achieving the specified level of performance under the constraints 
imposed pushes the frontier of current knowledge in control science and 
technology. In this section, we summarize our present understanding 
of—and contributions toward—a solution to the problems of active 
structural vibration suppression and precision optical axis regulation 
while accommodating wideband disturbances. The problem of flexible 
spacecraft reorientation is considered separately in Section 10. 





Because of the geometric complexity of the optical support 
structures being considered, a distributed-parameter model of the 
structural dynamics is generally not feasible. Although lumped-parameter 
models are readily obtainable via the finite-element method, quite high 
orders for such models (in terms of the number of vibrational modes 
retained) are needed to give an adequate representation of the structural 
dynamics. For the process of synthesizing an active controller, however, 
a reduction in model order is virtually unavoidable. Although the 
computational burden connected with carrying out the mathematical 
operations involved in a design is a consideration in this, the more 
fundamental reason is the unresolvable (deterministic) uncertainty in the 
computation of the higher characteristic structural frequencies. This 
exposes two matters with serious implications for control: 

(1) Parameters characterizing the full-order structural model are not 
known precisely; and (2) The reduced-order model used to represent the 
structure in the controller synthesis process Ignores potentially 
important dynamics, mostly at high frequencies. 


5.1.2 Initial Approaches to the Synthesis Problem 

For relatively small amplitude motions, the problems of structural 
vibration suppression and optical axis regulation are legitimately 
treated using linear dynamic models. Because of the geometric complexity 
of the structures, however, nonnegligible coupling of the structural 
modes through the actuators (and '■.he disturbances) is to be expected. A 
multi-input multi-output context for the controller synthesis is there¬ 
fore essential. Fortunately, the science of control for linear multi- 
variable systems—although still a fruitful research field—is quite 
mature [e.g., Refs. 5-1, 5-2, 5-3]. Many systematic design approaches 









Our initial approach to the synthesis problem consisted of 
applying several variants of the linear-quadratic-Gaussian (LQG) method 
for the determination of controller feedback strategy, each one assuming 
different degrees of availability of information about the disturbance. 
The object of control for the design exercises was (and remains) the 
generic ACOSS Model No. 2 (shown in Fig. 5-1) representing an optical 
support structure [cf. Refs. 5-4, 5-5]. Initial results, together with a 
description of the specific wideband disturbance affecting the structure, 
were given in Reference 5-6. With the feedback loop closed around a 
reduced-order (10 mode) model, each of the three initial designs 
exhibited stability; in addition, substantial reduction of the optical 
axis pointing error was affected (in fact, the performance specifications 
were satisfied when full knowledge of the disturbance statistics was 
assumed). These designs proved unstable, however, when the feedback loop 
was closed around an enlargement of the reduced-order design model 
incorporating 10 additional structural modes (mostly at higher 
frequencies). Without making unreasonable assumptions in regard to 
inherent structural damping, attempts to modify the basic LQG design so 
as to obtain closed-loop stability relative to the enlarged structural 
model (e.g., by incorporating an inner loop for velocity output feedback 
[Ref. 5-7], and/or by Incorporating actuator and sensor synthesizers 
[Ref. 5-8]) proved fruitless. We note that this sort of impasse is not 
an isolated occurrence [Ref. 5-9]; often, such Impasses are not highly 
publicized, however—perhaps unfortunately so. 

One of the principal conclusions we drew from this impasse was 
that not enough serious attention had been given to the questions of 
selecting reduced-order models or of selecting active transducers 
(actuators and sensors). In fact, these two processes are strongly 














Interdependent. Together, they have a decisive Influence upon the 
effectiveness of any controller feedback strategy. We also observed that 
whereas model selection, transducer selection, and (especially) feedback 
strategy determination have each received detailed study—in Isolation 
from one another—very little Is understood about the details of their 
Interaction. Among the reasons for this Is that such Interaction studies 
with models of significant complexity do not lend themselves well to 
theoretical inquiry. We therefore resolved to undertake an experimental 
approach (not in hardware) in the spirit of Reference 5-10 to Increase 
our understanding of these crucial interactions. This approach was 
presented in detail in Reference 5-11. In brief, it consists of a 
sequence of end-to-end synthesis trials—from model selection, through 
transducer selection, to feedback control design and evaluation—with 
one or more aspects of one of these principal factors in the synthesis 
process being changed at each trial. Individual trials aim not so much 
at the best possible design as at exposing some important feature of the 
interactions within the synthesis process. 

Progress with this experimental approach was reported in 
References 5-12 through 5-14. The emphasis there was placed primarily 
upon expanding our development of a particular approach to active 
transducer selection and upon exploring the influence of various 
transducer selections on the synthesis process—the controller feedback 
strategy being fixed (an LQG-strategy incorporating full knowledge of the 
disturbance statistics [Ref. 5-15]). Reduced-order controllers were 
obtained that exhibit stability in the presence of enlarged structural 
models containing modes well beyond the frequency extent of the reduced- 
order design model. 

In the remainder of this section, a final account of progress with 
the active control synthesis experiment is given. A detailed description 
of results obtained during the current reporting period is provided in 
Section 5.2. The principal contributions throughout the course of this 












contract toward the synthesis of vibration controllers for spaceborne 
optical systems are summarized and placed in perspective in Section 5.3. 
Finally, several research questions growing out of the work reported here 
that remain open and that have profound implications for future flexible 
space systems, including directed-energy systems, are discussed in 
Section 5.4. 

5.2 Current Results in Active Control Synthesis 

During the final reporting period under the contract, work toward 
active control synthesis for wideband disturbance accommodation focused 
on the following specific topics: (a) Modification of LQG compensation 
strategy to enhance stability robustness relative to unmodeled high- 
frequency dynamics; (b) Examination of the potential for independent 
(noncolocated) selection of sensors to improve sensing of disturbance 
effects; and (c) Incorporation of high-resolution spectrum estimation 
data-processing in the active control synthesis process to provide more 
accurate reduced-order modeling information. Primary emphasis was placed 
on the first of these three topics. 

5.2.1 Stability Robustness to Unmodeled High-Frequency Dynamics 


Initial study of interactions between principal elements of the 
synthesis process [Refs. 5-13 through 5-15] was done without changing the 
controller feedback strategy. In spite of this, a surprising degree of 
stability robustness to high-frequency modes excluded from the design 
model was observed, reflecting in part a compatibility between the 
selections of a reduced-order design model and of (colocated) active 
transducers. Stability in the presence of unmodeled modes up to a 
frequency of 14.8 Hz was demonstrated (compared to the design model 
bandwidth of 0.148 - 3.24 Hz). This fact generated optimism that a 
substantial improvement in stability robustness could be demonstrated by 
appropriate modification of the LQG disturbance-rejection strategy. A 
systematic approach to such modification (LQG/Loop Transfer Recovery) has 
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is the vector of outputs from the sensors; z = (c^, e L 0SY ’ e DEFOCUS 
is the vector of optical system line-of-sight (LOS) errors to be 
regulated; and B^: v*m, Bq: v*y, C P : txv, Cy: fcxv, D: 3xv are 
the geometric influence matrices associated with the actuators, 
disturbances, displacement sensors, rate sensors, and regulated 
variables, respectively. The superscript ”T" denotes matrix transpose 


The approach to actuator selection we have proposed [Refs. 5-18, 
5-12 (Part 1)] seeks to extract from a (relatively large) class of 
actuator candidates a (relatively small) subclass which retains 
sufficient control authority over the variables to be regulated to 
achieve the specified performance (with an appropriate feedback 
strategy). A sequence of least-squares approximations associated with 
equations having the form 


.LOSX 


(5-4) 


Ay 


(5-5) 


is computed, where the columns of A consist of actuator modal influence 
vectors 4> b^ : n *l from Eq. (5-1) representing candidate 
actuator locations, and the objective vectors are the (transposed) modal 
coefficients b L0SX = bL ° SY ” '^l.OSY 2 nx * tlle ^ ine “ 0 ^ -8i 8 bt 

error vector from Eq. (5-3). Successive iterations exhibit a monotonic 
decrease in the column dimension of A in Eqs. (5-4), (5-5) corresponding 
to a deletion of actuator candidates. (Criteria for decisions on the 
deletion of actuators is discussed in detail in Reference 5-12.) The 

iteration is terminated when an appropriate measure of either one of the 
LOSX LOSY 

least-squares residual vectors b - Ax, b - Ay exceeds a 
specified tolerance. One or more minimal selections (which do not exceed 
the tolerance parameter) are obtained. 







5.2.2.1 Algorithm for Sensor Selection 


In the light of this strategy for actuator selection, examination 
of Eqs. (5-1) through (5-3) suggests a completely analogous but 
complementary approach to the selection of sensors in which selections 
are keyed to a potential capacity for detecting disturbance effects. We 
seek to extract from a (relatively large) class of sensor candidates a 
(relatively small) subclass which retains sufficient information content 
about the effects of the disturbances upon the structure to provide for 
an appropriate level of disturbance rejection. To clarify this 
interpretation, we give a precise formulation of the sensor selection 
problem, referring to Reference 5-12 (Part 1) for details of its 
solution. The sensor selection may be defined using the following two 
sets of parameters: 

(1) A collection of sensor candidates, consisting of the modal 

T 1 1 

coefficient vectors 4> c J : nxl, where c J is the j-th row of 
the geometric sensor coefficient matrix Cp (for a displace¬ 
ment sensor) or C v (for a rate sensor) in Eq. (5-2); and 

(2) A collection i J T of disturbance objective vectors, consisting 
of the modal influence vectors 4> T b D ^ : nxl, where is the 
j-th column of the geometric disturbance influence 

matrix Bp in Eq. (5-1). 

To sharpen the succeeding discussion, the disturbance vector up is 
identified with the two-vector u^ = (u^y» u^)^ 

representing the two disturbance inputs denoted in Fig. 5-1, where the 
subscripts denote the specific nodes (nos. 37 and 46) of the 
finite-element structural model at which the disturbances are applied. 

The collection ,'j then consists of the two elements 

37 A T 37 46 A „T,46 

t = $ b D , t = $ b D 
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where bjj and are the columns of the matrix corresponding to the 
disturbance inputs U37 and u^g, respectively, in Eq. (5-1). The 
mathematical selection problem is a weighted-least-squares minimization 
corresponding to the algebraic equations 

Sx = t 37 (5-6) 

Sy - t 46 (5-7) 

where the columns of S are elements from the collection £P of sensor 
candidates to be considered. Several notational conventions will 
facilitate a precise statement of the sensor selection problem. Let X 
denote the number of candidates in a particular collection fP , and let 
•P' be a subcollection of SP . Denote by l(SP') the subset of (1,•••,!} 
that identifies which elements of , ( P are retained in SP' , and let S' be a 
matrix whose columns are the elements of SP '. 

Sensor Selection Problem . Find a subcollection SP' of SP and 
corresponding coefficient vectors x' = (xj), y' = (yj), j e l(SP') 
such that 

(1) J w (x'; S’, t 37 ) - (t 37 - S’x’) T W(t 37 - S’x’) and (5-8) 

J w (y’; S’, t 46 ) = (t 46 - S’y’) T w(t 46 - S’y’) (5-9) 

are minimized, where S-’ = \'P'\\ and 

(2) the subcollection SP' is minimal with respect to acceptable 
increases in the minimum values: 






(5-10) 



In the above, the weighting matrix, W: nxn is chosen to be real, 
symmetric, and positive definite. The problem statement is independent 
of the choice of column arrangement within S' [Ref. 5-12 (Part 1), 
Proposition 2-7]. 

Since the sensor selection problem has the same mathematical 
structure as the selection problem defined for actuators (cf. Eqs. (2-34) 
through (2-37) of Reference 5-12), the algorithm used to generate 
actuator selections may also be applied to sensor selection. Details of 
the algorithm have been given previously [Ref. 5-12, Sec. 2.2]. In 
brief, one begins by selecting an initial class of sensor candidates 
A sequence of least-squares minimizations of functionals having the form 
of Eqs. (5-8), (5-9) is then performed. At each stage of the sequence, 
the least-squares residual vectors just computed suggest one or more 
deletions of sensor candidates from for the succeeding stage 
according to guidelines indicated in the algorithm description. A 
lattice of eliminations is thereby generated. Termination of the 
sequence occurs when the minimum value in Eq. (5-10) or (5-11) 
(equivalently, an appropriate norm of the residual vector associated with 
either Eq. (5-6) or (5-7)) exceeds a specified tolerance. 


plications with ACOSS Model No. 2 


For a specified configuration of actuators, the length (however 
defined) of the least-squares residual vectors associated with 
Eqs. (5-4), (5-5) provides a measure of how well the actuators (as 


represented by the actuator modal influence vectors appearing as columns 
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and b L0SY ). In results reported previously [Refs. 5-13, 5-14], all of 
the actuator selections were generated so that this approximation was 
quite good (residual vector tolerance parameter = 10 -6 ). It should be 
noted in passing that the quality of this sort of static criterion of 
approximation can only be judged in the light of dynamic performance of a 
closed-loop design employing such actuators. 

However, in previously reported work, rate sensors were assumed to 
be colocated with the actuators selected without a corresponding (static) 
evaluation of the resulting sensor configuration. The sensor selection 
procedure presented in the preceding subsection provides a mechanism for 
making such evaluations. For any configuration of sensors, the length of 
the least-squares residual vectors associated with Eqs. (5-6), (5-7) and 
arising out of the minimization of the functionals in Eqs. (5-8), (5-9) 
provides a measure of how well the sensors (as represented by the sensor 
modal coefficient vectors appearing as columns of the S matrix) can be 
combined linearly to approximate the disturbances (as represented by the 
modal influence vectors t^7 and t^6). One may now ask how well (in 
the above sense) the (colocated) sensor configurations assumed in 
previous work approximate the disturbances. The answer—in some 
cases—is: not so good. Correlation of this observation with the 
performance and stability robustness of controller designs reported 
previously is interesting, as is indicated shortly. 

The specific assumptions under which design experiments have been 
conducted to date are briefly recalled. The structural model 
representing the object of control is Revision 1 of ACOSS Model No. 2 
[Ref. 5-5]. A modal damping ratio of Ci “ 0.001 is assumed for all 
modes (cf. Eq. (5-1). The model is subject to (wide-sense stationary) 
wideband disturbance excitation with frequency content as shown 
(one-sided) in Fig. 5-2. The disturbance is applied at two 

(uncorrelated) locations (cf. Fig. 5-1): node 37 on the upper optical 

support truss, and node 46 on the equipment section (Figs. 5-3 ff. 

display the node points explicitly). To hold down computation costs, it 





is further assumed that the points of application do not change with 
time, and that each disturbance is constant in direction (at an equal 
angle with each body axis). The specified optical performance is to 
regulate the optical axis (z-body-axis) rotation errors (LOSX, LOSY) 
about the two lateral axes (x and y body axes) so that 


LOSX _< 

0.05 V 

rad 

(rms) 

(5-12) 

LOSY < 

0.05 y 

rad 

(rms) 

(5-13) 


and to regulate the mirror relative positions so that a sharp image focus 
is maintained: 

DEFOCUS < 10" 3 m (rms) (5-14) 

where each inequality is in the root-mean-square (rms) sense. (Explicit 
expressions for the parameters LOSX, LOSY, and DEFOCUS as a linear 
function of nodal displacements are given in Reference 5-5.) 

Choice of reduced-order models for design is based upon ranking 
the response of the (uncoupled) modes of vibration for the open-loop 
(i.e., u A = 0 in Eq. (5-1)) structure to the total disturbance 
excitation. An abbreviated ranking is shown in Table 5-1. A corre¬ 
sponding interlaced reduced-order design model, together with a 
contiguous counterpart obtained by adding interleaved modes, is 
illustrated in Table 5-2. A family of variable-order contiguous models 
is used for evaluations of closed-loop stability and optical system 
performance with reduced-order controllers. 















Table 5-1. Ranking of open-loop modal responses to wideband disturbance 


RMS tOS-Error 
(u~rad) 


24 


600.97 
455.77 
23 o 105.89 

13 o 92.65 

22 o 77.04 

12 o 72.43 

16 o 32.47 

14 o 15.91 

10 o 12.94 

8 o 12.22 

5.06 


39 


5.01 


Rank Mode RMS l.OS-F.rror 
(u-rad) 


26 11 


42 20 


47 17 


78 18 


Table 5-2. Characteristic frequencies for AC0SS Model No. 2 and 


reduced-order models; o Design (interlaced) 

• Design (contiguous) ■ Evaluation (contiguous). 


Mode 



Frequency 

Mode 

Frequency 




01*) 


(list) 

7 

O 

• 

0. 148 

25 

■ 5.162 

8 

O 

• 

0.282 

26 

■ 5.170 

9 

O 

• 

0.319 

27 

■ 7.877 

10 

O 

• 

0.335 

28 

■ 7.918 

1 I 


• 

0.468 

29 

* .8.772 

12 

O 

• 

0.583 

30 

■ 8.776 

13 

O 

• 

0.601 

31 

■ 8.838 

14 

O 

• 

0.67 3 

32 

■ 8.980 

15 


• 

0.960 

33 

■ 9.620 

16 

o 

• 

1 .092 

34 

■ 10.40 

1 7 


• 

1.839 

35 

■ 11.67 

18 


• 

1 .344 

36 

■ 12.27 

19 


• 

1 .889 

37 

■ 13.39 

20 


• 

1 .990 

38 

■ 13.62 

71 


• 

2.060 

39 

■ 14.77 

7? 

o 

• 

2.452 

40 

■ 16.43 

71 

O 

• 

7.4 72 

41 

20.83 

24 

o 

• 

3.242 

4 2 

21.79 





4 3 

2 1.84 
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5.2.2.2.1 Evaluation of Sensor Configurations Defined by Colocation 

A summary of the minimal actuator selections generated in the 
context of the experimental approach to the synthesis process is 
presented in Table 5-3. The number of actuators p in a minimal selection 
is a function (among other things) of the number of modes n retained in 
the reduced-order design model, the functional type of the actuator 
candidates (.V) considered, and the tolerance parameter tq for the 
least-squares computation used to terminate the selection algorithm. The 
tolerance parameter is applied (as a maximum permissible upper bound) to 
the maximum (in absolute value) component of the current 

residual vector. The value of the tolerance parameter was allowed to 
fluctuate slightly (between 10“^ and 5 x 10“5) as an experimental 
variable. Note (Table 5-3) that for each selection: 


max LOSX 


max LOSY . 

i=l,...,n p i < T 0 


max DEFOCUS . 
i-l,...,n P i 1 T 0 


(5-15) 

(5-16) 

(5-17) 


although only Eqs. (5-15), (5-16) were enforced to terminate the 
selection algorithm. This is another reflection of the fact that 
Eq. (5-14) is the least stressful of the performance criteria Eqs. (5-12) 
through (5-14). The actuator locations for the minimal selections of 
examples A, B, C, E, and F of Table 5-3 are pictured in Figs. 5-3a, 5-4a, 
5-5a & 5-5b, 5-6a, and 5-7a, respectively (there being two minimal 
selections for example C). 

We now examine the sensor configurations that result from placing 
(rate) sensors of the corresponding type at the locations of the various 
actuator selections. Recall that each minimal actuator selection is 
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identified with a matrix A': nxp (where p is the number of actuators in 
the selection) and vectors x', y': pxl corresponding to least-squares 
solutions of Eqs. (5-4) and (5-5), respectively. The functionals 


V* 


defined by Eqs. (5-8) and (5-9), respectively, are then minimized, 
generating least-squares solutions x”, y”: pxl of Eqs. (5-6), (5-7), 
respectively, with the actuator matrix A' playing the role of a sensor 
matrix S. Maximum components of the corresponding residual vectors 
p 37 > p 46. n *i are listed in the last two columns of Table 5-3. The 
results show that least-squares solutions of Eq. (5-7) that are possible 
using the various colocated sensor configurations is generally poor, and 
most notably so using the minimal selection of example E (in which case 
the least-squares solution corresponding to Eq. (5-6) also represents a 
poor approximation). 

A correlation between these observations and dynamic performance 
under closed-loop control can be seen by examining Table 5-4. Data in 
this table summarizes design parameters for, and stability and 
performance obtained with, LQG/disturbance-rejection designs (not 
incorporating loop transfer recovery) that were done to examine the 
influence of reduced-order model selection and of active transducer 
selection upon the synthesis process. (Alphanumeric entries in the first 
column of Table 5-4 identify a spo. Lfic actuator selection from among the 
selections described in Tables 5-3, 5-5, and 5-6. For example, C-12a 
refers to the colocated 12-actuator selection a of Example C described in 
Table 5-3.) Performance data in Table 5-4 for example E using the 5 
colocated rotational actuators and (rate) sensors shows a marked 
degradation in the error parameter LOSY (and a lesser degradation in 








Controller design and evaluation; LOG/disturbance-rejec 
without loop transfer recovery. 

















LOSX) when a Kalman filter is connected to the LQ regulator in the 
feedback compensator. Since sensor measurements are inputs to the Kalman 
filter, one would expect less degradation upon incorporation of the 
filter through the deployment of sensors which provide better information 
on disturbance effects. This expectation is confirmed, as we demonstrate 
shortly. 


5.2.2.2.2 Independent Minimal Sensor Selections 

Using the algorithm for sensor selection described in Section 
5.2.2.1, several minimal selections for sensors have been generated. 

Data concerning them is summarized in Table 5-5, and they are Illustrated 
in Figs. 5-3b, 5-4b, 5-5c & 5-5d, 5-6b, and 5-7b (there being two minimal 
selections for example C). In generating these selections, the same 
reduced-order design model, the same functional type of transducer 
candidates, and the same tolerance parameter were used as for the corre¬ 
sponding actuator selections (cp. Table 5-3). Note that in example A, 
the number of sensors cannot be reduced below the number of modes in the 
design model without one of the residual vector components exceeding the 
tolerance parameter. Sensor selections with example C are rather 
remarkable. Comparison of Figs. 5-5b and 5-5c reveals that the locations 
and functional type of sensors in sensor selection a are identical to 
that of the actuators in actuator selection b. Stated differently, the 
sensor configuration defined by colocating (rate) sensors of 
corresponding type with actuator selection b is a minimal selection in 
its own right with respect to the disturbance environment. This 
remarkable fact is reflected in closed-loop LQG design performance. 
(Designs have only been done employing actuator selection a, together 
with colocated sensors, but reference to Figs. 5-5a and 5-5b and the data 
in Table 5-3 indicate only slight differences relative to selection b.) 
Table 5-4 indicates that an LOG controller (without loop transfer 
recovery) employing colocated actuator/sensor selection C-12a exhibits 
completely acceptable performance degradation associated with 
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incorporation of the Kalman filter into the feedback compensator, and 
displays the largest degree of robustness with respect to unmodeled high- 
frequency dynamics of any of the unrecovered designs for which p < n is 
assumed. 

Sensor selection with example E indicates a rapid deterioration in 
the least-squares approximations associated with Eqs. (5-6), (5-7) as 
sensors are reduced in number below n (the row rank of the sensor- 
candidate matrix S). In fact, the selection process terminates when p 
(the number of sensor candidates) is reduced to 15. This is consistent 
with the observations noted above regarding the poor least-squares 
approximation quality of the colocated sensors associated with example E 
(cf. Table 5-3). This suggests considering the (noncolocated) 
combination of the 5 actuators selected (Fig. 5-6a) with the 15 sensors 
selected (Fig. 5-6b). Note that this particular combination is 
noncolocated only in that the number of actuators and of sensors is not 
the same. As indicated in Table 5-4, the employment of this combination 
(designated E-5/15) in an LQG design leads to noticeable performance 
improvement relative to the colocated selection in example E. On the 
other hand, there is a slight degradation in robustness to unmodeled 
high-frequency dynamics associated with the additional sensors. The 
performance improvement (at the design model level) associated with the 
augmented transducer configuration is qualitatively predictable. The 
change in robustness to unmodeled dynamics is not. 

5.2.2.2.3 Further Observations on Actuator Selection 

The notation of a minimal selection (either of actuators or of 
sensors) is motivated in part by the potential for advantageous 
employment of nonminimal selections (to be obtained by the judicious 
enlargement of minimal selections). In particular, the process of 
convergence to a minimal selection gives sharp guidelines for such 
enlargement. Data for nonminimal selections referred to in Table 5-4 is 
given in Table 5-6. Design experience embodied in Table 5-4 indicates 
that the enlargement of minimal selections is not a predictable approach 
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to enhancing robustness to unmodeled high-frequency dynamics. Whereas 
the enlargement In example D of minimal selection D-9a to nonminlmal 
selection D-14 [cf. Ref. 5-12] results In enhanced robustness, just the 
opposite effect is observed in example E by enlarging minimal selection 
E-5 to nonminlmal selection E-lla. (Selection E-lla, illustrated in 
Fig. 5-8a, Is a heuristic augmentation of selection E-5 based on 
engineering intuition, and does not coincide with any of the several 
selections of 11 actuators—one of which Is depicted in Fig. 5-8b—that 
occur during the convergence of the selection algorithm to configuration 
E-5 (Fig. 5-6a)]. The role of enlargement to nonminlmal selections 
remains to be given systematic attention. 

The consequence of imposing appropriate constraints on the 
actuator (or sensor) candidate class needs to be understood at a deeper 
level for further exploitation of the selection process being applied 
here. It has been noted [Ref. 5-13] that selection E-5 (Fig. 5-6a) may 
not be feasible for an actual system since it implies an actuator at the 
center of mass of the primary mirror (node 1001). Example F (cf. Table 
5-3) gives the selection that results if the actuator candidate class 
excludes actuation devices from any of the mirror mass-centers. The 
number of actuators in the minimal selection is doubled (cf. Fig. 5-7a). 

A final remark is specific to ACOSS Model No. i,* The distinction 
between examples G and H (Table 5-3) reflects the subtl^fact that the 
degrees of freedom in the LOS-error model are not identical to the 
degrees of freedom associated with the kinematic mounts [cf. Ref. 5-5]. 
Specifically, the degrees of freedom of the secondary mirror are defined 
in terms of z-translation at, and x-axis and y-axis rotation about, its 
mass center, rather than in terms of the kinematically equivalent 
kinematic mount displacements. 







5.2.2.3 Remarks 


Although the strategy for Independent selection of sensors has 
sharpened insight regarding the effective deployment of active 
transducers, the full potential inherent in the selection process for 
selecting sensors to enhance irformation on disturbance effects has not 
been exploited. The initial step in the selection process executes a 
reduction in the number of candidates from X (the size of the candidate 
class to n (the number of modes retained in the reduced-order desigr 
model). The mechanism that realizes this reduction is the QR-factoriza- 
tion of the sensor candidate matrix S (cf. Eqs. (5-6), (5-7)), which 
generates a ranking among the sensor candidates based solely upon 
relative linear independence [cf. Ref. 5-12 (Part 1)J. Only in 
subsequent steps of the selection algorithm are eliminations of sensor 
candidates based upon least-squares approximation of disturbance 
influence vectors. 

On the other hand, there is a serious question as to the 
advisability of tailoring a selection of sensors too tightly to 
parameters of the disturbance environment which are not precisely known. 

A more fundamental course would be to investigate robustness in both 
stability and performance with respect to uncertainty or change in 
parameters of the disturbance environment as it relates to sensor 
selection. 

5.2.3 Incorporation of System Identification into the Synthesis Process 

As discussed in Section 5.1.1, there are fundamental limitations 
on the accuracy with which complex flexible structures can be represented 
using finite-element techniques. Moreover, there are no feasible 
alternatives to such modeling approaches at present. These facts provide 
strong motivation for seeking to enhance modeling capability through the 
incorporation of system identification techniques, especially for the 
specific purpose of increasing the accuracy of reduced-order models to be 
used for the synthesis of active controllers. The use of appropriate 










Identification techniques has the potential to reduce uncertainty in the 
knowledge of basic structural parameters (e.g., mass or stiffness 
distribution), and to provide a quantitative assessment of the input- 
output effect produced by delays and small-amplitude nonlinearities 
associated with the dynamics of active transducers, the sampling and 
reconstruction of analog signals, and digital computation. 

The high-resolution spectrum estimation procedure originally 
formulated by Fogel [Ref. 5-19] and extended by Villalba [Ref. 5-20 and 
Sec. 8 herein] is very attractive for application to reduced-order 
modeling of flexible spacecraft for several reasons. First, the 
mathematical formulation accommodates the essential constraints on data 
availability consistent with a flight environment. It assumes a 
capability for system excitation at the location of each actuator and at 
the location of each anticipated disturbance, and the availability of 
measurement data at each sensor location. Second, the procedure is aimed 
specifically at solving the problem of identifying modes that are 
closely-spaced in frequency—an inherent property of large flexible 
structures. In this subsection, we describe the data transformations 
required in order to incorporate this identification approach into the 
synthesis process, and point out a fundamental limitation which requires 
some adjustment of the synthesis process itself. 


5.2.3.1 Data Transformation Required for Incorporation 


The high-resolution spectrum estimation procedure is itself a 
dynamic process that observes the input-output behavior of the system to 
be identified, and generates the spectral characteristics (natural 
frequencies and damping ratios) of a linear time-invariant system which 
exhibits the same input-output behavior. It is assumed that actuator and 
sensors have been selected, both by functional type and by location, and 
that the locations of expected disturbance inputs to the system are 
known. Specific input data required by the identification process 
consists of discrete-time samples of all possible time histories of the 









y r (t; u k ) , 0<t£t 1 ; k - 1,... ,m , r - 1, (5-18) 

y r (t; w k ) , 0 < t < tj ; k-l,*M,Tf (5-19) 

where y r (•; u^) represents the response function at sensor location 
r to an impulsive excitation at actuator location k, and y r (»; w^) 
represents the response function at sensor location r to an impulsive 
excitation at disturbance location k (cf, Eqs. (5-1), (5-2)). 

Specific output data provided by the identification process are the 
numerical parameters in the modal representation for the estimate of the 
responses observed: 


y r (t! V ■ l" te [ (a kr + J b ki> Plt ] 

i-1 


»i - - Vi + i »i (1 - <i 2 >‘ 


n is the estimate of the number of vibration modes, and u>^ and are the 
estimates for the natural frequency and damping ratio of the i-th 
estimated mode. In what follows, the letter j represents the imaginary 
unit /-I. 

In order to use the information contained in Eqs. (5-20), (5-21) 
for reduced-order modeling of the observed system, it is necessary to 
generate related estimates for the matrix products associated with the 
eigenvalues (i.e., 0^ an d 2Z0) and eigenvectors (i.e., $ T B^, * T Bp, Cp$, 
Cyi, and D$) that appear in the dynamic equations Eqs. (5-1) through 
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(5-3). (Note that explicit generation of the eigenvalue matrix ft and the 
eigenvector matrix 4 themselves is not needed.) It is immediately clear 
from elementary algebra that the entries of ft 2 = diag(o>i 2 ) and 
2Zft = diag (2$^) are given, respectively, by 

“ Re 2 (p t ) + Im 2 (p 1 ) 

and 

2 ?i“l ■ “ 2 Re(p ± ) 

where Re(») and Im(*) represent the real and imaginary parts, 
respectively, of the indicated arguments. Generation of the eigenvector 
matrix products is somewhat more involved and leads to the need for 
imposing certain mild but subtle restrictions upon the synthesis process. 


5.2.3.2 Deduced Estimates for Eigenvector Matrix Products 

Estimates for the eigenvector matrix products are obtained by 
comparing the expressions Eqs. (5-20) to the analytic representations for 
the various Impulse responses, which are easily written. First, it is 
quite helpful to develop a concise yet descriptive notation for 
representing the entries of the various matrix products. In the finite- 
element context, the matrix 4: v*n has the form 

[4 1 : ••• : 4 n ] (5-22) 

where the i-th column vector consists of components <J>^ that 
represent the amplitude of the i-th mode shape at node o of the finite 
element model and in the kinematic degree of freedom B (e.g., rotation 
about one of the body axes) at that node. The number of such distinct 











pairs (a,8) represented in the model is the total number v of generalized 
coordinates. The geometric actuator influence matrix B A : v*ra has the 
form 


where each column vector represents a particular actuator, with a 
structure dependent upon the actuator type. For an actuator whose effect 
can be fully described at a single node (e.g., rotation about one of the 
body axes at that node), the corresponding vector b^ is a unit vector of 


with the 1 at the index (a,8) defining the relevant node and degree of 
freedom. For an actuator involving a connection between two nodes, the 
corresponding vector b£ is a unit direction-cosine vector with non-zero 
entries at the appropriate pairs (o',6'), of indices. 

Consequently, the modal actuator influence matrix $ T B A : nxm can be 
represented in the form 


if each actuator represents nodal actuation or rotation. If one colum 
(say the r-th) represents an axial node-connecting actuator, then the 
r-th column has the slightly different structure 






v''CvIXV.wf VLVl -.1 
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C r) 


where each entry is a sum 

♦(act r) ’ *!♦(»',p') " ♦(o-.B-)l co ‘ - e '• “"• 6 " ) 

involving the direction cosines of the node-connector. The subscript 
"(act r)" denotes the index (a r ,8 r ) defining the node and degree of 
freedom corresponding to the r-th actuator. The notation of Eq. (5-23) 
is used repeatedly in what follows. A corresponding representation for 
the modal sensor influence matrices Cp4> and Cy$ may be given; e.g.. 


The modal disturbance influence matrix $ Bp is represented similarly 
to Eq. (5-23) with the entries being denoted ^ instead. The 

structure of the modal line-of-sight coefficient matrix Cp 0S $ ■ D$: 3*n 

has no correspondingly simple representation. The entries are in general 

. . . - . .i _LOSX .i „LOSY . .i .DEFOCUS 

inner products of the form <p • C , <f> »C , and $ *C , 

respectively, involving the full mode-shape vectors if* 1 (cf. Eq. (5-22)) 

and the coefficient vectors in the LOS-error model [Ref. 5-5]. 

Analytic expressions for the various impulse responses with the 
system of Eqs. (5-1), (5-2) are presented next, employing the notation of 
Eqs. (5-23), (5-24). The details of the development are elementary and 
are therefore omitted. It is assumed that the zero-frequency modes are 
not represented (i.e., * 0) and that the system modal damping ratios 

Ci satisfy 



0 < <1 


( 5 - 25 : 


To avoid obscuring the principal ideas, results are presented only for 
the employment of rate sensors (i.e., Cp ■ 0 in Eq. (5-2)) and for the 
employment of nodal actuators and sensors (i.e., Eqs. (5-23), (5-24) 
apply without change). The general case follows directly (though 
tediously). Under assumption Eq. (5-25), the r-th (rate) sensor response 
to a unit impulse at the k-th actuator is given by: 


7 T <*1 u k > " I ♦?, 
r K 1-1 


(5-26) 


V t; u k ) 

and 


y (act k) e H u, i t 
“d. 


[(-?lWl)sin(u> d t) + u) d cos(Uj t)] (5-27) 


(5-28) 


After some algebra, Eqs. (5-26) through (5-28) may be placed in the form 
of Eq. (5-20) (with the carets on and n removed). The 

parameters a* r » b* r produced by the identification process are then seen 
to have the decomposition: 


\r 


9 (v-sen r) T (act k) 


♦fv-i 


(v-sen r) 9 (act k) - 1/2 

(1 - cf) 


(5-30: 
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Noting that Eq. (5-30) is a scalar multiple of Eq. (5-29), the problem of 
constructing the desired matrix estimates reduces to solving Eq. (5-29) 
for the two (sets of) factors. 


5.2.3.2.1 Estimates for the Actuator and Sensor Influence Matrices 

We rewrite Eq. (5-29) emphasizing the associated array of indices: 




(v-sen r) 9 (act k) 


(5-29a) 


corresponding to the totality of impulse responses. It appears that 
Eqs. (5-29a) do not point to a unique solution for the individual factors 
unless some restriction is imposed on the degree of noncolocation of the 
actuators and sensors. (For example, a single actuator and sensor— 
noncolocated—lead to a one-parameter family of solutions). Uniqueness 
of solutions can be demonstrated, however, under the following mild 
restrictions: 


(Al) The actuator and sensor configuration contains at least one 
colocated actuator/sensor pair (both by node and degree of 
freedom); and 

(A2) For each mode i - l,»»»,n, there is a colocated actuator/ 

sensor pair k, for which the parameter a, , in Eq. (5-20) 

ii 

for the response estimate y. (•; u, ) is nonzero. 

“i i 

Three cases are distinguished. 


Case 1: Colocated actuator/sensor pairs . The appropriate subset of 
Eqs. (5-29a) corresponding to colocated actuator/sensor pairs is 


i 

a kr 


v (act r) v (act k) 


1,•••,M , r - 1,...,M ;!-!,• 


(5-29b) 


90 










where M < min {m,i} is the number of colocated pairs. For each mode 
1 - l,...,n, solution for the factors may be obtained from the diagonal 
subsystem of Eqs. (5-29b): 


*(act k) 9 (act k) 


, k - !,••• ,M 


and Is easily seen to be 


4 141 (5 - 3I) 


where the caret "■»” denotes a (deduced) estimate of the Indicated 


Case 2 : Noncolocated (rate) sensors . The appropriate subset of 
Eqs. (5-29a) corresponding to sensors which are not colocated with any 


1 A 1 » 1 

a kr “ ^(v-sen r) *(act k) 


, k - !,•••,M , r - M+l,•••,£. ; 1 - 1»••• »n 

(5-29c) 


Fix a mode 1 in {l,«**,n}. Estimates $* act , k - 1,***,M, have been 

determined in Case 1 by Eqs. (5-31). By assumption (A2), there is a 

number k. In for which af . * 0. It follows from Eqs. (5-29c) 

1 x i K 1 

that <^ act ) * 0, and therefore that 


*(act k ) 


K^J 172 (, k 1 k 1 ) 


, r - ttfl,•••,*. 








Case 3 : Noncolocated actuators . The appropriate subset of Eq. (5-29a) 
corresponding to actuators which are not colocated with any sensor is.: 


'•’(v-sen r) '•’(act k) 


, k = , 


• ,M ; i = 1, • • • ,n 
(5-29d) 


Fix a mode i in {l,»»*,n>. The sensor factors <t l ( v _ sen r ) appearing in 
Eqs. (5-29d) represent sensors colocated with actuators, for which 
estimates have been determined in Case 1 by the identification: 


together with Eqs. (5-31). By assumption (A2), there is a number r^ in 
{1,•••,M} for which a* #0. It follows from Eqs. (5-29d) that 

Vi 

_ x* 0, and therefore that 


I 1 i 

a sgn (a ) 

I Vi 


Taken together, Eqs. (5-30) through (5-32) constitute the desired 
estimates for the matrices $ B A and C v $ in Eqs. (5-23) and 
(5-24), respectively. 


5.2.3.2.2 Est imates for the Disturbance and LOS-Error Influence Matrices 

The analytic representation for the disturbance response is 
analogous to Eqs. (5-26) through (5-28). Under assumption Eq. (5-25), 
the r-th (rate) sensor response to a unit impulse at the k-th disturbance 
location is given by: 
















reduced-order model. It can be assumed therefore that at least one 
component, say $( v _ gen r )» °^ the vector 1° Eq. (5-38) is nonzero. Then 
from Eq. (5-36), 1 


^(v-sen r ) 


(5-39) 


where the appropriate explicit expression (either Eq. (5-31) or (5-32) is 

to be used for <j>^ . in Eq. (5-37). This establishes an estimate 

(v-sen r .) 
r 1 


for the matrix $ Bp. 

An attempt to estimate the LOS-error Influence matrix 
(cf. Eq. (5-3)) leads to a certain difficulty which directly affects the 
synthesis process. Recall that this matrix has the structure: 


<w> 


(^.C L0SX ) 



(* n .C L0SX ) 

(♦".C L0SY ) 



In order to estimate the entries In C^os*! estimates of the modal 


vector components <J> 


(o,8) 


at each index (o,6) which is represented In the 


LOS-error model are required. However, estimates of such components are 
available only at indices (o,8) corresponding to locations of actuators, 
sensors, or disturbance inputs (cf. Eqs. (5-31), (5-32), (5-33), (5-39)). 
In general, there is no guarantee that any of the required components 


helps to clarify this fundamental point. Consider the colocated 
actuator/sensor selection E-5 (cf. Tables 5-3, 5-4) consisting of 5 
transducer pairs. The Indices (o,8) corresponding to actuator (and 
sensor) locations are (cf. Figs. 5-6a, 5-6b): 
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where 8(.) denotes the rotational degree of freedom about the Indicated 
body axis. Indices corresponding to disturbance input locations are 


(o,8) - (37, x), (37, y), (37, z), (46, x), (46, y), (46, z) 

where the second element in each pair represents the translational degree 
of freedom along the indicated body axis. It can be verified, however, 
[Ref. 5-5], that none of these index pairs appear in the LOS-error 
model. Therefore, for this example, the best estimate of the matrix 
Clos $ that can be computed (cf. Eq. (5-40)) with the available data is 
zero! 

5.2.3.3 Implications for the Synthesis Process 

The inability to generate feasible estimates for C LO s* bears 
upon the active controller synthesis process as follows. In the 
LOG/disturbance-rejection designs used to generate controller feedback 
strategy, the functional being minimized is [cf. Ref. 5-15]: 

J ■ / (z T Qz + p R uj u A ) dt (5-41) 

where z: 3x1 is the variable to be regulated (LOS-error vector) appearing 
In Eq. (5-3). Therefore the matrix product Clqs* is a factor (twice) 
in the product z T Qz. The use of poor (or zero!) estimates for 
c L0S $ in E< 1* (5-41) generally leads to poorly performing controllers, 
since no penalty for large LOS-errors is imposed In the design process. 

There are several alternatives available to resolve the present 
impasse to the Incorporation of the high resolution spectrum estimation 
technique into the synthesis process. Th^ functional Eq. (5-41) may be 



altered so as to weight only the sensor measurements of the system 
Eq. (5-1) (i.e., z is replaced by the sensor output y in Eq. (5-41)). 

More directly, the LOS-error vector could be treated as a sensor output 
Itself, an apparently reasonable assumption for certain optical systems 
[Ref. 5-21]. Experience with LQ controller designs in which the 
LOS-error vector is not weighted in the cost functional Eq. (5-41) 
indicates difficulty in meeting the performance specifications [Ref. 5-6; 
Section 5], Therefore, the latter alternative appears preferable. 


Since neither of these alternatives has been implemented in the 


synthesis experiments performed to date, a database for evaluating the 
effect of incorporating the particular identification process under 
discussion into the synthesis process is not available. However, the 


process itself has been exercised on ACOSS Model No. 2 yielding excellent 


Principal Contributions 


The principal contributions toward an understanding of the active 
control synthesis process that are embodied in References 5-11 through 
5-15 and 5-18, together with Sections 5 and 6 of the present report, are 
briefly summarized. 


Several generally understood principles have been given a sharp 
focus in the specific context of wideband disturbance accommodation with 


precision flexible space structures. First, the inherent inseparability 
of reduced-order model selection, active transducer selection, and feed¬ 
back strategy determination within the total synthesis process has been 
clearly demonstrated. In particular, the consequences of inadequate 
attention to the first two factors cannot be overcome by even the most 


sophisticated feedback strategy. Second, the value of • 


approach to the study of the mutual interactions of these three factors 
has been established. In particular, a thorough understanding of the 
synthesis process is equally as important as (and the key to generation 
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of) a specific design product. Third, the fundamental compromise between 
stability and performance (more precisely, between stability robustness 
and performance robustness) with respect to unmodeled high-frequency 
dynamics has been clearly exhibited. In particular, the results obtained 
highlight the need for more systematic design tools to ensure performance 
robustness for high-precision flexible structures. 

In addition, specific contributions have been made toward 
techniques of active transducer selection and controller feedback 
strategy determination. 

Active Transducer Selection . First, a new and systematic approach to the 
selection of actuators and sensors has been developed. Selection of 
actuators is related directly to parameters of specified system 
performance, and a distinct selection of sensors may (if desired) be 
related directly to parameters of the anticipated disturbance 
environment. Selection algorithms are based on theoretically sound 
principles; reliable portable software for implementation is widely 
available. Second, the effectiveness of the selection strategy in the 
active control synthesis context has been demonstrated. In particular, 
the role of active transducer selection within the overall synthesis 
context has been clarified. 

Controller Feedback Strategy Determination . First, the effectiveness of 
a systematic approach [Refs. 5-16, 5-17] for enhancing stability 
robustness to unmodeled high-frequency dynamics has been convincingly 
demonstrated in a stressful and realistic design context relevant to 
flexible structure control. Second, important limitations of this 
approach in obtaining performance robustness with high-precision flexible 
structures have been revealed. 
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Having given an account of the work accomplished, it is 


appropriate to discuss significant matters that remain unresolved. There 


are two categories: matters relating specifically to the work 
accomplished; and matters of general concern. 


We begin with some specific matters. First, the strategy for 


active transducer selection that has been developed has been shown to 


work surprisingly well. However, a theoretical explanation for this 
success is lacking. Are the transducer selections optimal in some 


well-defined sense for the purposes of actuation or of measurement? 


are the "static" criteria for selections linked with the closed-loop 


dynamic performance associated with the resultant selections? Answers to 


these and related questions would provide a key to much more efficient 
implementation of the selection strategy. For example, the strategy 
seems to be potentially valuable for predicting the consequences upon 
performance of prior constraints on transducer selection, or of 
transducer failures. Second, the emphasis to date in feedback strategy 
determination has been upon robustness to unmodeled high-frequency 
dynamics. Although we believe this to be a proper priority, robustness 
to model parameter variations is the logical—and no less important— 


next step. Third, the way has been cleared for the incorporation of 


high-resolution spectrum estimation into the reduced-order modeling 


process. Significant investigations should include the development of a 
data base of design experience to assess the observable effects upon 


stability and performance of such incorporation, and to assess the effect 
of accounting for additional delays in the signal path (e.g., due to 


actuator and sensor dynamics, or to data sampling, digital computation, 
and data reconstruction). 


Finally, we touch upon several matters of overriding general 
concern. First is the matter of robustness in performance. Most 


systematic approaches to robustness enhancement for the feedback control 










of systems with (deterministic) uncertainty focus principally upon 
robustness of stability. Guidelines for augmenting performance 
robustness typically amount to suggesting a tightening of the original 
performance specification. For very high precision systems such as the 
present application, such guidelines are not feasible. The need for 
effective guidelines to assist in augmenting performance robustness lies 
in the following fact: if a system exhibits (deterministic) uncertainty 
relative to its design model representation, performance is likely to 
deteriorate more sharply than stability. From an operational viewpoint, 
a stable system that is unable to accomplish the assigned mission is—at 
best—a waste of resources. The point is that fundamental research on 
this topic is needed. Second is the matter of robustness (in both 
stability and performance) with respect to disturbances. Disturbances 
that can be anticipated for real systems often defy precise 
characterization, both as to their analytic structure and parametric 
description, as well as to the location and direction of their influence 
upon the system. The consequences of inadequate characterization of 
disturbances, as well as techniques for enhancement of robustness with 
respect to disturbances are highly significant, and unresolved, 
questions. 
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SECTION 6 


STABILITY ENHANCEMENT IN STRUCTURAL CONTROL 


6.1 Introduction 

In our experimental approach to active control design for flexible 
space structures [Ref. 6-1], numerous actuator and sensor placements have 
been explored [Refs. 6-2, 6-3]. The main goal of the accompanying 
controller designs [Ref. 6-2] for these transducer selections has been to 
ensure that the state-feedback portion of each linear-quadratic-Gaussian 
(LOG) design met the line-of-sight (LOS) performance specifications as 
closely as possible. For each corresponding LQG design the largest set 
of modes, contiguous in frequency, that produced a stable evaluation was 
noted. In most instances to date closed-loop stability was maintained 
with modes well beyond the 18-mode design model, but stability with all 
150 elastic modes of the structural model was not achieved. 

This section presents sufficient conditions in the frequency 
domain for improved stability robustness to unmodelled high frequency 
dynamics and outlines a design procedure to enhance the robustness of our 
present LQG designs. We demonstrate these ideas on three of our 
previously reported controller designs, and indeed achieve closed-loop 
stability for two of them in the presence of all 150 elastic modes of the 
structural model. We conclude with observations on our controller 
results in the context of the experiment, and discuss some of the issues 
for flexible space structure control that have yet to be addressed. 









6.2 Present Controller Results 

The control design methodology used with the various 
actuator/sensor selections in connection with our experimental approach 
Is detailed in Reference 6-2. In brief, the methodology involves the 
synthesis of an LQG controller for a design model consisting of the first 
18 elastic modes (modes 7-24) of Revision 1 of ACOSS Model No. 2, 
augmented with a model of the two wideband disturbances which Impact the 
structure at the points indicated in Fig. 5-1. A statistical description 
and state-space representation for the disturbances can be found in 
Reference 6-2. In the performance index for the linear-quadratic (LQ) 
regulator design the three LOS outputs, LOSX, LOSY (optical axis rotation 
errors) and defocus, are weighted with a diagonal matrix to ensure that 
the state feedback controller meets the LOS specifications [Ref. 6-2] as 
closely as possible. A Kalman filter is then used to estimate the states 
from the sensor outputs. Table 6-1 lists controller design results 
reported on thus far. Included in the table are LOS output weightings, 
'll* <1 2 > ^3» and the scalar control weighting, p R , used in the LQ 
performance index. The scalar parameter, pg, which multiplies an 
identity matrix to represent the sensor noise covariance and ensure a 
nonsingular Kalman filter problem, Is also listed. Relevant performance 
and stability results are noted In the table. 

Up to this point in the experiment, the main focus has been the 
achievement of adequate LQ-regulator performance with different 
actuator/sensor selections, while the stability of the LOG controllers 
has been an observed variable. Even though no particular effort has been 
made to enhance whatever stability robustness properties our standard LQG 
designs [Ref. 6-2] may possess, the number of modes in the stable 
evaluations extends well beyond that of the 18-mode design model. 

However, closed-loop stability has not been achieved with all 150 elastic 
modes of the structural model, so the robustification of our designs must 
be of prime concern. 
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Fortunately, the lack of stability robustness has a clearly 
identifiable cause. Specifically, the only perturbation in our design 
model that we employ to test the robustness of our controllers is the 
addition of 132 modes, all of whose frequencies are above those of the 18 
modes in the design model. Other perturbations, for example, variations 
in structural parameters which could change the frequencies and mode 
shapes of all 150 modes, have not yet been treated. Robustness to 
parameter variations is an important issue itself. The goal we pursue 
here is to modify the LQG controller design to enhance its robustness to 
high frequency modes not in the design model, preferably without using 
explicit knowledge of those modes. 


Frequency Domain Conditions for Stability Robustness 


In this subsection we present conditions which guarantee 
closed-loop stability of a control system when it is evaluated with high 


frequency dynamics excluded from a design model of the plant. We first 
discuss frequency domain descriptions of transfer matrices, and then 
describe appropriate matrix expressions for design model errors, as back¬ 
ground to deriving the criteria for stability robustness. These 
frequency domain criteria motivate the design modification to our 
controller, presented in Section 6.4, intended to enhance stability 


robustness. 


In our investigation of robustness to high frequency dynamics, we 
examine closely the frequency domain properties of the transfer matrix 
HCsJPjCs). Here H(s) represents the transfer matrix for the 
compensator and P<j(s) is the transfer matrix for the 18-mode design 
model of the plant (i.e., the structure to be controlled) without 
disturbance dynamics appended. The transfer matrices are shown in 
Fig. 6-1 as part of a closed-loop multi-input multi-output (MIMO) control 
system. The output y, the control u, and the reference variable r are 








The product H(8)Pd(s), which is sometimes referred to as a 
return ratio matrix, or a loop transfer matrix, is the matrix multiplying 
signals inserted in the control loop broken at the input to the plant, 
point x, to produce the signals at the output of the same breaking 
point. Similarly the loop can be broken at the output of the plant, 
point xx, to produce a corresponding loop transfer matrix P<j(s)H(s). 

For single-input single-output (SISO) control systems, no difference 
exists between the transfer functions P^H and HP<j (with the Laplace 
variable s suppressed in an obvious way) because they are both scalars. 

In general, for MIMO systems the products corresponding to different 
breaking points will be different because they are transfer matrices 
(matrix multiplication being generally noncommutative). Our reason for 
focusing on the loop transfer matrix HP<j corresponding to breaking the 
loop at the input to the plant will become apparent later. 

The frequency response properties of a SISO system are determined 
from familiar gain and phase plots. For MIMO systems, analogous 
frequency domain information is available about the gain through the use 
of matrix norms. Here we make use of the matrix spectral norm. For any 
matrix A, the maximum singular value, o(A), and the minimum singular 
value, o(A), are defined as: 
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°< A) “ ^min 1 ^ 1 ) 


where X[•] denotes an eigenvalue of the indicated matrix and the super¬ 
script H denotes complex conjugate transpose. 

It is well known that these two singular values have the following 
important characterization: 


a(A) 



( 6 - 1 ) 


£(A) 


min 

x*0 



( 6 - 2 ) 


Equations (6-1) and (6-2) indicate that a and £ represent, respectively, 
the maximum and minimum input-output gain ratio, viewing A as a transfer 
matrix [Ref. 6-4J. Therefore maximum and minimum singular values may be 
taken as appropriate measures of frequency response for a transfer 
matrix. 

Before we develop our robustness criteria, we also need an 
appropriate matrix description of the error in the design model of the 
plant as compared to the evaluation model. As we have indicated, the 
error we are concerned about is the exclusion of high frequency modes 
from the design model. The most straightforward way of describing that 
error follows from partitioning the transfer matrix for the full order 
model (in our case the evaluation model): 
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P f (s) 


[°d °h] 


(si - A d ) _1 0 


B d 

0 (si - A h ) _1 _ 


- B h - 


■ v st -v‘ B d + c h ( 8 i -v'\ 

- P d (s) + P h (s) (6-3) 

where subscript "f" refers to the full order model, subscript "d” refers 
to the design model, and subscript "h" refers to the neglected high 
frequency modes. 

Equation (6-3) is an additive description of the model error and 
is presented here to introduce the concept and notation associated with 
such error representations. An alternative manner of describing this 
same error yields the most useful information on how the loop shape HP d 
should behave in the frequency domain for stability robustness. In 
particular, we use the so-called multiplicative form of the error, 
reflected at the input to the plant and written as: 

P f (s) - P d (s)[l + L(s)] (6-4) 

where 

L(s) & P’^s) [P f (s) - P d (s)] - P d ! (s) [P h (s)] 

As indicated in Reference 6-5, the frequency domain behavior of 
L(s) as measured by the variation with frequency of its maximum singular 
value, a(L), is smaller than unity in the low frequency region, but much 
larger than unity at high frequency where Ph(s) has significant 
dynamics. 

With this background material we can now proceed to determine 
sufficient conditions for closed-loop stability when evaluating our 
compensator on a full-order model. 
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Because the focus of our concern is the control loop broken at the 
input to the plant, we first examine the closed-loop transfer matrix 
between the reference variable r and the control variable u at the plant 
input of Fig. 6-1: 


■ H(s)P.(s)j H(s) r(s) 


The stability of the closed-loop system can be examined via the 
determinant det[I + HP<j]. In particular, if 4 >ol( s ) represents the 
characteristic polynomial of the open-loop transfer matrix HP<j and 
<|>Cl( 8 ) represents the characteristic polynomial of the closed-loop 
system, then [Ref. 6-7]: 


* (s) det[l + HP.] 


From this relationship the multivariable Nyquist criterion can be 
derived [Refs. 6-6, 6-7]. This criterion says that the closed-loop 
system of Fig. 6-1 is stable if and only if the number of clockwise 
encirclements of the origin, 0, by the image of the clockwise Nyquist 
contour, D R , under the mapping det [I + HP<j] plus the number of 
unstable poles of the open loop transfer matrix, HPj, is equal to 
zero. In symbols, this is sometimes written: 


N(0, det[I + HP d ], D r ) + Py ( 


By designing the compensator, H(s), via the LOG methodology we are 


guaranteed that the closed-loop control system of Fig. 6-1, with 18 
design modes of the plant, will be stable, which means of course that the 


multivariable Nyquist criterion is satisfied. The LQG methodology does 
not guarantee, however, that a closed-loop system in which the 














compensator is connected with a plant other than the design model will be 
stable. Using the multivariable Nyquist criterion, we determine 
sufficient conditions which guarantee that closed-loop stability is 
maintained when the full-order plant, Pf(s), is connected with the 
compensator; in symbols: 

N(0, det[l + HP f ], D r ) + P U>QL * 0 (6-5) 

where here the open loop transfer matrix is HPf. 

Our discussion in part follows the presentations of References 6-8 
and 6-6 and that of Reference 6-5. Consistent with Theorem 2.2 of 
Reference 6-6 and Theorem 2 of Reference 6-8 we assume that: 

(a) The number of unstable poles (including those on the ju axis) 
of HP^ is the same as the number of unstable poles of HP^; 

(b) All poles of HPf on the jw axis are also poles of HP^; 

(c) The closed-loop system designed for P<j is stable. 


By employing the LOG design methodology and by using design and 
evaluation models that include only elastic modes of the space structure, 
we are guaranteed that conditions (a) through (c) are satisfied. 


With these assumptions we can focus on the conditions which 
guarantee that the encirclement count of the map of det[I + HPj] does 
not change as the determinant is perturbed to det[I + HPf], In fact, 
if the determinant is never zero as it goes from det[I + HPj] to 
det[I + HPf], then the encirclement count of the origin will not 
change, and Eq. (6-5) will hold. 


Using our matrix error description of Eq. (6-4), we can 
parameterize the deformation of the resulting Nyquist diagram so that 
closed-loop stability is ensured. The condition is that 


det[l + HP (I + (L)] * 0 


(6-6) 







for all € in the interval [0,1] and all s on the contour Dr. 

Because both maximum singular values and a(P^) approach zero 

as | s j ■+■ <*>, we only need to verify Eq. (6-6) for values of s along the 
portion of the contour D R lying in the closed left-half plane, i.e.: 

J2 r = {s: s e D r and Re(s) 0} 

Presuming that the inverse of HP^ exists, we can rewrite 
Eq. (6-6) as: 

det[(HP d )(l + (HP d ) -1 + C L)] * 0 
The condition to be satisfied then reduces to: 

det[l + (HP d ) _1 + C L] * 0 (6-7) 

for all € on the interval [0,1] and all s on the contour 

Equation (6-7) is not considered a good measure of how close the 
perturbed system, with error represented by L(s), is to being unstable, 
because the determinant of a matrix is a poor indicator of how close that 
matrix is to being singular. A more reliable measure of nearness to 
singularity is provided by matrix singular values. With this in mind we 
note that a general matrix sum A + E is guaranteed to be nonsingular, or 
equivalently to have det[A + E] * 0, if 

o(A) > o(E) 










so the condition of Eq. (6-7) will be true if 

a[l + (HP d ) _1 ] > e a(L) 

for all e on the interval [0,1]; and in particular if 

o[l + (HP d ) _1 ] > o(L) (6-8) 

Equation (6-8) is a well-known sufficient condition for stability derived 
from the multivariable Nyquist criterion which is valid if the design and 
evaluation models are related by Eq. (6-4). 

Using the minimum singular value inequality: 

o(I + A) < o(A) + 1 

we have for Eq. (6-8): 

! + £[ Hp d ] -1 > oJl+OlP^" 1 ] > o(L) (6-9) 


for all 0 < u < -. 

In our earlier discussion of multiplicative model error we noted 

that 

o(L) » 1 









, equivalently: 


-— > o(L) 

“[ Hp d ] 


5[hp ] < 


It follows that 


o[HP d ](u)) « 1 (6—1C 

for those frequencies such that [a(L)](u>) » 1. Equation (6-10) 
establishes the design goal that is used to judge the robustness of oui 
control system designs. It is the MIMO generalization of the SISO 
requirements of low loop gains in the frequency range of unmodelled 
dynamics. 

6.5 Design Procedure to Enhance Stability Robustness 

Now that we have identified the cause of instability in our LOG 
designs and have established a design goal to help alleviate that 
problem, we discuss in this subsection how our LOG controllers can be 
modified to more closely satisfy Eq. (6-10). 

To illustrate the problem at hand explicitly, examine Fig. 6-3a 
which is a plot of the maximum and minimum singular values (cr-plot) of 
the loop transfer matrix HP,j for actuator/sensor selection E-5 of 
Table 6-1, consisting of five rotational actuators and sensors. At no 
frequency less than 100 Hz (mode 88 is at 99 Hz) does the maximum 
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singular value cross below 0 db. Such frequency domain behavior is 
certainly not in the spirit of the design goal represented by Eq. (6-10), 
and indeed Table 6-1 indicates that the system does become unstable when 
mode 34 (10.4 Hz) is added to form an evaluation model of modes 7-34. 

The a-plot of Fig. 6-3a is for the transfer matrix of our LQG 
design where estimates of the state are being fed back to the 
linear-quadratic regulator gain G. Consider now the situation where the 
actual states are being fed back, namely that of the LQ regulator which 
is depicted in closed-loop form in Fig. 6-2. 



Figure 6-2. LQ regulator. 


By breaking the loop at the plant input we get the open-loop transfer 
matrix for an LQ regulator which is G$(s)B, where $(s) denotes the 
resolvent matrix (si - A)~l associated with the state-variable 
representation of the structural dynamics: 

x = Ax + Bu 

The a-plot associated with the L0 loop for the state-feedback design 
employing 5 rotational actuators is shown ,in Fig. 6-3b. Clearly this 










state feedback open-loop transfer matrix, with a crossover at 30 Hz, has 
frequency domain properties more in line with our design goal than does 
the LQG loop matrix HP d . Our estimator, as reflected in HP d , has 
ruined what was a more desirable loop shape when the actual states were 
being fed back. This degradation in loop shape caused by using a Kalman 
filter to estimate the states is apparently quite common [Ref. 6-9] and 
is observed in all of our examples. 

Recent studies [Refs. 6-9, 6-5] have shown that the loop shape of 
the LQG loop matrix HP d can be made to approach the loop shape of its 
corresponding LQ regulator G$(s)B by the use of a full-state loop 
transfer recovery procedure. Such a procedure is desirable of course 
because our LQ loop shapes have lower gains at high frequency than those 
of HP d . 

The derivation of the recovery procedure can be found in 
References 6-9 and 6-5. Here we only sketch the details. 

We first require the design model of the plant, which in our case 
is augmented with a disturbance model, to have outputs greater than or 
equal in number to the inputs. The design model should also be minimum 
phase, though the recovery procedure is possible with a nonminimum phase 
plant (except that recovery will generally be poor at frequencies 
corresponding to the nonminimum phase zeros). 

To achieve full-state loop transfer recovery, the Kalman filter is 
redesigned using a modified algebraic Riccati equation for the 
determination of the filter gain: 


where the subscript "a" on the state, control, and observation matrices 
A, B, and C, respectively, makes reference to the fact that the design 
model of the plant is augmented by a model of the disturbance [Refs. 
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Figures 6-3. 


a-plots for designs with 5 rotational transducers 
elastic modes 7 through 24. 
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6-10, 6-11], and the matrix B a(1) is the Influence matrix of the white 

noise source upon the augmented state. The new term In Eq. (6-11), 

T 

qB a B a , can be Interpreted as representing an additive fictitious 
process noise. By allowing q -*■ <*>, It can be shown [Ref. 6-9] that 
asymptotically, the Kalman gain, K, satisfies: 

* B W(p I) _1/2 (6-12) 

/q 

where W is an orthogonal matrix. It is shown in Reference 6-5 that with 
the proper substitution of Eq. (6-12) into the expression for the open 
loop transfer matrix (where the subscript "a" is suppressed for simpli¬ 
fication) 

H(s)P d (s) - (G[sl - A - BG - KC] -1 K)C$(s)B 

along with appropriate matrix algebra, that indeed: 

HP., G$B 
a 

polntwise in s, though not uniformly, as q increases to infinity. 

Our discussion up to this point has concerned loop shapes 
corresponding to breaking the loop at the input to the plant (cf. 

Fig. 6-2). The same robustness results can be stated for loops broken at 
the output of the plant, where in this case the relevant loop transfer 
matrix is P^H. Likewise the full-state loop transfer recovery procedure 
has a dual for loops broken at the plant output. For such a recovery 
procedure, the LQ regulator is modified so that the open loop transfer 
matrix, P^H, approaches the Kalman filter loop transfer matrix, C$K. 

This is a desirable course to follow presuming that the Kalman filter 
loop has frequency domain properties superior to those of P<jH. To 









explore this alternative a little further, examine Figs. 6-3c and 6-3d 
which show the Kalman filter loop shape and the loop shape of PjH, 
respectively, for the controller designs employing five rotational 
actuators and sensors. In this case P^H has smaller loop gains in the 
high frequency region than does the Kalman filter loop, so recovery of 
the filter loop would be counterproductive. 

Because our LQ design has been specially tailored to ensure 
adequate LOS performance, and because the LQ loop already provides 
desirable high frequency behavior, we have chosen, for all the examples 
to be discussed, to modify only the Kalman filter via the full-state loop 
recovery procedure so that the LQG loop matrix HPj approaches the LQ 
loop matrix G$B asymptotically. In the traditional jargon associated 
with the active control of flexible space structures, our focus on 
achieving low gains at high frequency for loops broken at the plant input 
is similar to minimizing the "control spillover". 


Demonstration of Stability Enhancement 


This subsection Illustrates how the full-state loop recovery 


procedure may be used in achieving our goal of low loop gains at high 
frequency for compensators connected with three distinct actuator/sensor 


placements on ACOSS Model No. 2. Some of the practical problems 


associated with the entire design procedure are noted. 


The desirability of using full-state loop transfer recovery for 
the compensator associated with the five rotational actuators and sensors 
has been discussed in the previous subsection. A value of q * 1 * 10^ 
for the recovery procedure was found to give numerical problems. With 
q = 1 x 10^ in the adjustment procedure, excellent stability robustness 
results were achieved. First examine Fig. 6-3e which is the a-plot of 
HPfj after the recovery procedure. This is virtually an exact copy of 
the target LQ loop shape of Fig. 6-3b, indicating good full-state 
recovery. Indeed an evaluation of this new compensator on all 150 




elastic modes of the structural model produced a stable closed-loop 
system. So the compensator, which had previously maintained closed-loop 
stability in the presence of modes 7-33, has been adjusted to maintain 
stability with modes 7-156. 

For this new LQG compensator design, closed-loop performance when 
the compensator is connected to the 18 design modes was tested (see 
Table 6-2). Compared to the LQG/18-mode performance results of the 
original compensator, the errors in LOSX and LOSY increase substantially. 
Reference 6-9 clearly Indicates that values of q sufficiently high result 
in filters that do a poor job of noise rejection, and this is manifesting 
itself in the poorer performance results we have noted here. Because the 
original L0G/18-mode controller (q *» 0) also violates the performance 
specification, a later design example was deemed more appropriate for 
investigating potential performance improvements by lowering the q value. 

Our design procedure is next demonstrated on the control system 
originally designed to employ nine translational actuators and sensors, 
denoted D-9a in Table 6-1, and originally discussed in Ref. 6-2. The 
frequency domain behavior of the LQG loop matrix HPj and the LQ loop 
matrix G$B for this design is shown in Figs. 6-4a and 6-4b respectively. 
The maximum singular value <3 of the LQ loop crosses the 0 db line at 
35 Hz while the maximum singular value a of the LQG loop HPj crosses 
the 0 db line well beyond 100 Hz. Apparently the LQ loop presents 
desirable target dynamics for recovery and hence for stability robustness 
enhancement. A value for the recovery parameter, q, of lx 10 ^ gave 
numerical problems, but values of 1 x 10? and lower gave acceptable 
results. Figure 6-4c shows the recovered LQG loop matrix, HP<j, 
indicating good recovery and a 0 db crossover at 35 Hz. 

However, when this new compensator is evaluated on modes 7-62 of 
the structural model, five pairs of poles occur in the right-half plane. 
These unstable poles include a pole corresponding to mode 30 and poles at 
frequencies between those of modes 34 and 38. These results indicate 
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Figure 6-4d. LQ loop transfer matrix G$B: redesigned for 
stability robustness, o R = 10”^. 
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Figure 6-4e. LQ loop transfer matrix GIB; redesigned for 
LOS performance, modified Q. 























that this design is particularly sensitive to the addition of high 
frequency dynamics, and even the crossover at 35 Hz provided by the LQ 
loop is not sufficient to bring any improvement to the original 
compensator. 

To explore how this stability problem might be alleviated, the LQ 
regulator was redesigned employing a much larger control weighting, 

Pg = 1 x 10~13. Reference 6-5 indicates that by raising the control 
weighting, the entire LQ loop shape tends to be lowered in magnitude, 
which Implies that the 0 db crossover would occur at a much lower 
frequency. This indeed is the effect as indicated by Fig. 6-4d, which 
shows 0 db crossover occurring at 6.1 Hz for the new regulator design. 
Note that this a-plot is shown between only 1 Hz and 100 Hz. The 
compensator, H(s), is completed by designing a Kalman filter using 
Pq * 1 * 10”*3 and full-state loop recovery with q = 1 x 10?. An 
evaluation on modes 7-62 produces closed-loop stability. 

The price that is paid for this improved stability result is a 
decrease in performarice, which was checked for the state-feedback portion 
of the design. This poorer performance is indicated by the higher 
root-mean-square (RMS) LOS error shown in Table 6-3 for the LQ regulator 
using pr = 1 x 10“13 compared to the errors produced by the regulator 
using pr - 1 x 10 ” *7 also shown in Table 6-3 for actuator/sensor 
selection D-9a. 

Because acceptable LQ LOS-performance is considered a minimum 
requirement, we need to modify this regulator design. To achieve LOS 
performance within the specification, the Q matrix weighting the LOS 
outputs in the performance index is modified to: 










while retaining the control weighting * 1 x 10 . The performance of 

the regulator is thereby brought within the specification and is shown in 
Table 6-3. The a-plot of the LQ loop for this modified design is shown 
in Fig. 6-4e and indicates a 0 db crossover at 35 Hz, exactly the same 
crossover frequency of the original design. 

To complete the compensator design, the Kalman filter with 
Pq = 1 * 10“13 is modified via the full-state recovery procedure 
using q = 1 x 10?. An evaluation with modes 7-62 produces the same 
five pairs of unstable poles noted earlier. These results indicate an 
apparent conflict between the objectives of meeting the LOS performance 
specifications with the LQ design and meeting stability robustness 
requirements of the LQG design for this selection of nine translational 
actuators and sensors. 

Such stability/performance conflicts are common in control system 
design. One could continue iterating on the compensator design, but the 
approach more likely to bring success is to examine a new actuator/sensor 
selection. This is in the spirit of our "experimental" approach. 

The final controller design examined is that employing twelve 
translational actuators and sensors (selection C-12a in Table 6-1). The 
LQ and LQG designs (prior to the use of loop transfer recovery) for this 
transducer selection achieve performance and stability results quite 
favorable in comparison to most of the other designs (see Table 6-1). 

The striking characteristic of the a-plot for the LQG loop matrix 
HPj shown in Fig. 6-5a is that the maximum singular value remains at a 
virtually constant value, well above the 0 db line, for most of the high 
frequencies of the plot. The a-plot for the LQ loop matrix shown in 
Fig. 6-5b indicates crossover at 51 Hz. 

An appropriate value of q for the full-state loop recovery is 
1000; values higher than this result in numerical problems. The 
frequency domain response of HPj is shown in Fig. 6-5c after applica- 









tlon of the adjustment procedure. When the new LQG compensator is 
evaluated on all 150 elastic modes in the model, closed-loop stability is 
maintained. 

For this new LOG compensator design, performance is tested when 
the compensator is connected to the 18 design modes (see Table 6-4). 
Compared to the LQG/18-mode performance results of the original 
compensator, the RMS LOS error has increased substantially and even 
violates the performance specification. As noted earlier, sufficiently 
high values of q in the recovery procedure result in Increased LOS 
error. We expect, therefore, that lowering the value of q should result 
in better performance, but may decrease our stability robustness. 

Indeed, a value of q *= 0.1 in the adjustment procedure leads to excellent 
performance results for the "recovered" LQG compensator evaluated with 
the 18-mode design model (see Table 6-4). Closed-loop stability with the 
full 150-mode model is also maintained. 

While we have demonstrated acceptable performance at the design 
model level and stability robustness, we must also examine performance 
robustness, i.e., whether adequate performance is maintained in the face 
of high frequency dynamics. An evaluation of the compensator (q ■ 0.1) 
with modes 7-36 produces unacceptable degradation in the LOS errors as 
indicated in Table 6-4. Not reflected in the table is the fact that the 
major contribution to LOS error comes from the disturbance located at 
node 37 on the upper portion of the optical support structure. 

We believe that a remedy for the poor performance robustness must 
begin with a reexamination of the modes in the design model, somehow 
taking into account those modes not presently in the model that neverthe¬ 
less make a substantial contribution to open-loop LOS error. Careful 
selection of actuators and sensors, and then careful design of a 
compensator for stability and performance robustness must follow. 
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Figures 6-5. o-plots for designs with 12 translational transducei 
elastic modes 7 through 24. 
















Recovered LQG loop transfer matrix HP 






















6.6 Discussion, Conclusions and Recommendations 

In the foregoing we first described how our control system designs 
for the numerous actuator/sensor placements on the space structure model 
are geared almost entirely toward achieving acceptable performance with 
the state feedback (or LQ) portion of our compensators. The observed 
lack of stability robustness of the LQG designs is the result of high 
frequency dynamics not included in the design model of the plant. 

Second, we derived from the multivariable Nyqulst criterion a design goal 
of low loop gains at high frequency as a means of enhancing stability 
robustness. Third, we discussed a design technique, full-state loop 
transfer recovery, that is a systematic aid to achieving the stated 
design goal, given that the LQ regulator has desirable high frequency 
properties. Finally, we demonstrated our revised design procedure on 
three compensators and noted the problems that can arise in practice. In 
particular, we saw for the design with nine translational transducers 
that the frequency domain requirements for LQG stability robustness 
conflict with the achievement of adequate performance even at the design 
model level of the LQ regulator. For the other two compensators, the 
recovery procedure allows us to achieve lower loop gains at high 
frequency and even stability with all 150 elastic modes. In every 
design, excessive values of the recovery parameter q produce numerical 
problems. The upper bound for q that can be tolerated must be determined 
for each design. We also noted in the final design that a value of q 
could be adequate for enhancing stability robustness, but may be 
sufficiently high that performance is severely degraded even at the 
LQG/design-model level. Again q will have to be adjusted for each design 
to achieve the desired results. 

We have not achieved performance robustness with our designs and 
have suggested that a remedy must begin with a reexamination of the 
design model of the plant. 
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In the broader context of our experiment we can see that reduced- 
order model selection, active transducer selection, and compensator 
design must not be considered in isolation. On the one hand, no matter 
how carefully a design model and active transducers have been selected, a 
poor controller design may preclude even the achievement of stability. 

On the other hand, a robust control system design cannot make up for an 
inadequate design model, or ill-considered actuator and sensor selections 

One of the major robustness questions we have yet to address 
explicitly is robustness to variation of parameters in the model of the 
structure. In the frequency domain context we require high loop gains in 
the frequency band containing the design model — usually residing in the 
low end of the frequency spectrum — to achieve robustness to parameter 
variations. Other significant questions for control design include the 
addition of actuator and sensor dynamics in our models and the 
examination of digital controllers. 

Clearly the stringent performance requirements in the face of 
significant disturbances acting on a flexible space structure continue to 
pose a formidable challenge in active control synthesis. 
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SECTION 7 


MODELING FOR FLEXIBLE STRUCTURE CONTROL WITH 
ELECTROMECHANICAL ACTUATOR DYNAMICS 


7.1 Introduction 

There are a variety of possible actuation devices that could be 
selected for use in achieving vibration control of large flexible space 
structures. Our recent research [Refs. 7-1, 7-2] has focused on electro¬ 
mechanical devices as actuators for achieving such vibration control. In 
the present section, we clarify our approach to the incorporation of 
electromechanical actuator dynamics by placing it in a unified context, 
and distinguish its principal features from other relevant work. 

An electromechanical actuator (also called a voice coil actuator) 
consists of two mass elements and an associated electrical circuit. As 
current is caused to flow through the circuit due to an applied actuator 
input voltage, an electromagnetic force is generated between the two 
actuator mass elements. A schematic is shown in Fig. 7-1. 



Figure 7-1. Basic electromechanical actuator. 
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Our motivation for use of electromechanical actuators Is as 
follows. Electromechanical actuators are simple, reliable, and adaptable 
to use as actuators for a wide variety of distributed structures in many 
different configurations. They can be used as linear (i.e., translation¬ 
al) actuators to generate control forces; and they can be used as rotary 
actuators to generate control torques. 

7.2 A Generic Case: Models for Possible Actuator Configurations 

Consider the control of the axial motion of a bar with distributed 
mass and distributed stiffness. There are three conceptually different 
structure/actuator configurations. These three configurations are 
presented and the general form of mathematical models used to describe 
the dynamics in each case is given. The derivation of the models follows 
a Lagrangian approach similar to that presented in References 7—1/ and 7-2 
and is therefore omitted. To clarify the main issues we wish to 
emphasize, it is assumed that inherent structural damping is negligible 
and that mechanical motion of the actuator does not of itself (i.e., 
kinematically) give rise to significant forces on the actuator. The form 
of the models is considerably more general in that they arise in cases 
where the structure is more complicated than a simple bar. Our 
consideration of this simple physical case should give an indication of 
the fundamental differences in the three configurations. In each 
configuration, certain simplifications of the mathematical models are 
indicated; in each case, conclusions are stated in terms of the overall 
way in which the actuator dynamics affect the structural dynamics. 

7.2.1 Structure and Inertial-Based Actuator 

Assumptions for this case are as follows. One actuator element is 
fixed in an inertial frame; the other actuator element is connected 
directly to (or is part of) the structure. A schematic of the generic 
case is shown in Fig. 7-2. 
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Figure 7-2. Representation of an inertial-based actuator. 

The dynamic model of the structure, for this configuration, can be shown 
to be given by 


Mq + Kq = Bf a 

where q is the vector of generalized structural displacements, M and K 
are generalized mass and stiffness matrices, respectively, for the 
structure, and B is the influence matrix of the actuator forces upon the 
structure. The electromagnetic actuator force f a on the structure is 
linearly related to the circuit current I: 

f =?I (7-1 

a 

which in turn satisfies the electromechanical coupling equation: 


a 


LI + RI (■ CB q 


(7-2 











where L and R are the inductance and resistance parameters, respectively, 
of the circuit, and u is the actuator input voltage. 


There are a number of reasonable simplifications of this model. 

For example, the inductance parameter L is typically small; if L is 
assumed to be zero, then the model can be described by 

Mq + 5 2 BR _1 B T q + Kq = £BR _1 u 

Thus, under appropriate assumptions, the dynamic characteristics of the 
structure are modified by the coupling with the actuator dynamics. In 
the special case that L = 0 is assumed, a structural damping effect is 
introduced into the structure dynamics through the actuator "back emf". 

Since there is no inertial base in a space environment, this 
configuration would not be a possible implementation to achieve vibration 
control of a flexible structure in space. The configuration is mentioned 
here primarily for comparison with the other two configurations to be 
discussed next. 


7.2.2 Structure and Structure-based Transmission-type Actuator 

Assumptions for this case are as follows. Each actuator element 
is connected to the structure through an elastic "transmission 
device". A schematic of the generic case is shown in Fig. 7-3. 

The dynamic model of the structure, for this configuration, can be 
shown to be given by 

Mq + Kq = BK (z - B T q) 
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Figure 7-3. Representation of a structure-based transmission-type actua 

and the model of the actuator elements is given by 
M a* z * + K ,< z - ' B a f a 

where z is the vector of generalized displacements of the actuator, M a 
and K a are generalized mass and stiffness matrices, respectively, for 
the actuator, and B a is the influence matrix of the actuator forces 
upon the actuator elements. The electromagnetic actuator force f a 
satisfies an equation of the form Eq. (7-1), where in this case the 
circuit current I satisfies 

LI + RI + 5bJz = u (7-3 

instead of Eq. (7-2). 


These equations are inherently coupled through the actuator flexi 
billty term. The effective actuator force on the structure is given by 








which can be viewed as a force transmitted from the actuator mass to the 
structure. It is this view that motivates the designation of the 
actuator as being of transmission type. Note that the above equations 
for the structure and the actuator are intrinsically coupled; it is not 
possible to view the actuator dynamics as a pure series connection with 
the structural dynamics. 

As with the preceding configuration, we look for reasonable 
simplifications of the model. The inductance parameter L is typically 
small; if L is assumed to be zero, then the model can be described by 

Mq + Kq = BK fl (z - B T q) 

M a *z‘ + £ 2 B a R _1 B a Z + K a (z - B T q) = $B R^u 

If, in addition, the stiffness of the "transmission device" is 
sufficiently large so' that it is reasonable to assume that K a + ®, then 
the model can be described by 

(M + BM a B T )q + ? 2 BB a R _1 B a B T 4 + Kq - £BB R -1 u 

Thus, under appropriate assumptions, the dynamic characteristics of the 
structure are modified by the coupling with the actuator dynamics. Even 
in the limiting case where L - 0 and K a + <*> are assumed, there is a 
modification of the structural mass characteristics and a structural 
damping effect is introduced into the structure dynamics through the 
actuator "back emf". 


7.2.3 Structure and Structure-based Reaction-type Actuator 


Assumptions for this case are as foilows. One actuator element is 
connected directly to (or is part of) the structure; the other ac 
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element is free to act as a "proof mass". The two actuator elements are 
connected through a "reaction spring". A schematic of the generic case 
is shown in Fig. 7-4. 
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Figure 7-4. Representation of a structure-based reaction-type actuator. 


The dynamic model of the structure, for this configuration, can be shown 
to be given by 

Mq + Kq - -BM fl (z + B T q) 

and the model of the proof-mass actuator element is given by 

M (z + B T q) + K z - B f 
a a a a 

where (in contrast to the preceding configuration) z is the vector of 
generalized relative displacements between the two actuator elements, and 
M a» K a» and B a are associated with the proof-mass element of the 
actuator. The electromagnetic actuator force f a satisfies Eqs. (7-1) 
and (7-3) (with the appropriate revision in the interpretation of the 
z-variable). 








These equations are Inherently coupled through the actuator 
Inertial term. The effective actuator force on the structure is given by 

f = -M a (*z + B T q) 

which can be viewed as a force in reaction to the acceleration of the 
proof-mass element of the actuator. It is this view that motivates the 
designation of the actuator as being of reaction type. Note that the 
above equations for the structure and the actuator are intrinsically 
coupled; it is not possible to view the actuator dynamics as a pure 
series connection with the structural dynamics. 

Again we look for reasonable simplifications of the model. The 
inductance parameter L is typically small; if L is assumed to be zero, 
then the model can be described by 

Mq + Kq = -BM a (z + B T *q*) 

M (z + B T q) + 5 2 B R _1 B T z + K z - ?B R _1 u 

a a a a a 

If, in addition, the mass of the "proof mass” is sufficiently large so 
that It is reasonable to assume that M a + », then the model can be 
described by 

Mq + S 2 BB a R -l B a B T q + (K + BK a B T )q = -^BB^^u 

Thus, under appropriate assumptions, the dynamic characteristics of the 
structure are modified by the coupling with the actuator dynamics. Even 
in the limiting case where L = 0 and M a + <® are assumed, there is a 
modification of the structural stiffness characteristics and a structural 










damping effect is introduced into the structure dynamics through the 
actuator "back emf". 

7.3 Feedback Control Issues 

By explicitly considering the structure and actuators together as 
the system to be controlled, our viewpoint is that the control input is 
the voltages applied to the actuator input circuits. 

There are many possible feedback control laws that could be 
examined as a means of achieving some (vibration) control objective. 

Only a few cases are mentioned here. 

a. The simplest control law is to assume that u = 0. In such a 
case the actuators are essentially used as passive dampers. 

b. Suppose that u depends only upon feedback of z and z; then 
only internal feedback of the motion of the actuator elements 
is required. In such a case the actuators are essentially 
used as "active" passive dampers. 

c. The general case is to assume that u depends on feedback of 
both the motion of the actuator elements and the motion of the 
structure. 

Certain aspects of the general case have been considered in 
References 7-1 through 7-4. The emphasis in References 7-1 and 7-2 (and 
in the present section) is on the development of general models for 
transmission-type and reaction-type devices as control actuators for a 
flexible structure. In the case where the structure is controlled by 
multiple actuators, a special form of decentralized control law is 
suggested and analyzed. A special form of centralized control law is 
also briefly examined. 

In contrast, the emphasis in References 7-3 and 7-4 is on the 
practical development of a linear reaction-type actuator. No attention 
is given to the development of a general model which characterizes the 
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coupling between the structure and the actuator. Control Issues are 
considered solely In the context of the developed device. 

Although there is a considerable literature on the design of 
control laws for achieving vibration control of flexible structures, the 
developed control laws must be suitably modified If the dynamics of 
electromechanical actuators are to be taken into account. 


7.4 Summary 

If electromechanical devices are used as control actuators for a 
distributed flexible structure, it is claimed that the dynamics of the 
actuator should be carefully examined as an Integral part of the 
structural system to be controlled. Apparently, there are several 
different structure/actuator configurations which differ in their 
physical characteristics and in the form of suitable mathematical models 
for these configurations. 

The important feature, from a control design standpoint, Is that 
it is relatively easy to Incorporate any actuator dynamics Into the 
structural dynamics In such a way that existing control design 
techniques, such as those considered In Sections 5 and 6 of this report, 
remain valid. 
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SECTION 8 


POLE LOCATION MANIPULATION FOR 
IMPROVED SPECTRAL RESOLUTION 


8.1 Introduction 

The problem of obtaining hlgh-resolution power spectrum estimates 
has been of concern to researchers In the area of signal processing for a 
very long time. This is due to the wide range of technical problems with 
solutions that are made feasible by the availability of sut_L estimates. 

A classic application of high resolution techniques is in the area of 
radar Doppler processing. Another application which happens to be of 
current relevance is the system identification problem for large flexible 
space structures. Because such structures characteristically exhibit 
vibrational modes which are very close in frequency, any identification 
procedure based on power spectrum analysis must employ a spectrum 
estimation technique with the ability to resolve closely spaced spectral 
peaks. 

In Reference 8-1, the problem of improving the resolution of power 
spectrum estimates was discussed. In particular, an estimation procedure 
was considered which is based on parametric analysis of the autocorrela¬ 
tion sequence {r(k)> corresponding to the wide-sense stationary discrete¬ 
time process (x(k)} of interest. The procedure uses the estimates 
(r(d), r(d+l), •••, r(d+N-l)} of the actual autocorrelation sequence 
(hereafter we shall assume d = 0) to obtain estimates of the parameters 

{A,, •••, A } and fz,, •••, z } of the following assumed structure for 
1 1 n '1 n ' 









It was demonstrated In Reference 8-1 that the locations of the poles 
{zj, •••, z n | in the complex plane have a significant effect on the 
sensitivity of the parameter estimates to errors in the autocorrelation 
estimates {r(0), r(l), •••, r(N-l)} (see Guideline 1, p. 211, and 
Guideline 2, p. 223 in Reference 8-1). Since the accuracy of the 
parameter estimates determines the accuracy and therefore the resolution 
of the resulting spectrum estimate, the implication is that if pole 
locations could be manipulated, then resolution could be Improved. 

In Reference 8-1, the use of desampling as a means for manipu¬ 
lating pole locations was considered, and a modified power spectrum 
estimation procedure based on desampling was constructed. In this 
section the improved resolution obtainable by the desampling procedure 
is demonstrated. But first, we consider the use of the autocorrelation 
sampling period as another means of manipulatir ; pole locations when 
discrete-time techniques are to be applied to the spectrum estimation 
problem for continuous-time processes. 


Sampling Consideratir 


for Continuous-Time S 


Consider a real, wide-sense stationary, continuous-time p 
and let its autocorrelation function r c> which is defined by 

r (t) - Efx (t) x (t + t)] 


epresentable by the analog of Eq. (8-1) for continuous-tim 








Here {C^, •••, C^} and {pp •••, p^} occur in complex conjugate pairs, 
and due to physical constraints we have Re(p^) _< 0, i = l,***,n. Since 
the power spectral density function Rj, of the process x c is 


? c i p i 

1 . 1 (jn - PjXJi! + p t ) 


a spectrum estimate may be obtained by estimating the parameters 

fCj, •••, C n f and {p^, •••, p R } (assuming n is known). This may be done 

by applying the spectrum estimation procedure from Section 7.2 of 

Reference 8-1 to the sequence (r(k)} obtained by periodically sampling 

r C : 


r(k) = r c (kT) , — < k(integer) < +» (8-4) 

where T is the autocorrelation function sampling period. Combining 
Eq. (8-4) with Eq. (8-2) we get 

n p Tiki 

r(k) = ? C e 1 1 1 , -a. < k < +» (8-5) 

i=l 


Comparing Eq. (8--5) with Eq. (8-1), note that 


P i T 

A t = C t , z t = e , i = 1,•••,n (8-6) 


Letting 


o t = Re(p t ) , ^ = Im^) , i = !,•••,n 


and assuming T is chosen such that 







-IT < fi.T < 


i - 1, 


(8-7) 


we may calculate {A^, •••, A n } and {z^, 
estimate of the Indicated variable) via 
and then obtain estimates {C^, •••, C n } 


• • • ,n 

• • •, z } (where M *“ denotes an 
n ‘ 

the procedure of Reference 8-1, 
and {pj, •••, p n } as follows: 


i = 1, • ••,n 


where Arg z^ denotes the principal argument (i.e., the unique real number 
which satisfies z = |z| exp(j-Arg z) and -it < Arg z £ it for any nonzero 

complex number z) of z^. Note that if Arg z^ » it, then in light of 

Eq. (8-7), z^ should be considered a numerical artifact and ignored by 
forcing - A^ - 0. This may be done by eliminating the columns of 
and V^, and the elements of and corresponding to A t and 
before solving Eq. (7-18) of Reference 8-1. 

Referring to Eq. (8-6), we see that there is an explicit relation¬ 
ship between z^ and T. Clearly, T may be used to manipulate the 
position of z^ in the complex plane. Furthermore, if 0 £ » 0, 
i - l,«»«,n, then z^ * i * l,«»«,n, and Guidelines 1 and 2 

may be used to judge the effect of T on the sensitivity of the spectrum 
analysis procedure to errors in the estimates (r(0), •••, r(N-l)}. If 
n is even and the set {$!.,•••, } is ordered such that 

0 < fijT < fi i+1 T < tt, i - 1, •••, n/2, we may set 

u 1 = « i T , i - 1, •••, n/2 
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(see Eq. (7-30) in Reference 8-1) and then use Guidelines 1 and 2 to 
determine values of T which will decrease the sensitivity of the spectrum 
estimation algorithm, implying an improvement in the resolution of the 
resulting spectrum estimates. 

In this subsection, emphasis is placed on the use of Guideline 2 
to improve resolution. Guideline 2 states that it is desirable to have 
uj c ■ tt/ 2 , where 

+ "n/2> I st 1 + a n/2 

“c 2 ' \ 2 

Since estimates for (Hj + &n/2)/2 are readily available in a number 
of applications (for example, radar Doppler processing), frequently T may 
be chosen on the basis of prior knowledge such that u> c = tt/2 . Thus 
Guideline 2 is a practical tool for improving resolution. 

The use of Guideline 2 is demonstrated here via three examples. 

In all three cases, the power spectrum of a continuous-time process x c 
is estimated via the algorithm of Section 7.2 of Reference 8-1. The 
required estimates for the samples r(k) = r c (kT), k * 0, •••, N - 1, 
are obtained from the data sequence 

(y(k) = x(k) + w(k), k = 0, •••, N<j - l} 
where x(k) denotes the sample 

x(k) ** x^(kT) 



of the continuous-time process x c , and w(k) is a Gaussian, zero mean, 
white noise sequence with variance o w 2. The autocorrelation 
estimator used is 


r(k) 


T y( O y( n+k) 




H 


( 8 - 8 ) 










which presumes the ergodicity of the process x c . Note that since 


E[yU) yU + k)] 


i r<k) + % 

^ r(k) , 


and because in this section it has been assumed that d ! 


equation in Eq. (7-16) of Reference 8-1, and the first and (N + l)st 
equations in Eq. (7-18) of Reference 8-1 were deleted prior to obtaining 
{a., •••, a } (the autoregressive coefficient estimates) and 


{A , •••, A } for the examples which follow. 


For the first example we have N<} * 128, N * 30, o w / = 0.1, and 


cos(0.9xt) + cos(irt) 


Nine cases are considered by choosing T so that u c = (0.95ir)T varies 
from 0. In to 0.9ir in increments of O.lir. For each case, ten power 
spectrum estimates are calculated, with each corresponding to a different 
realization of the noise sequence {w(k)}. Figures 8-la through 8-li 
contain plots of these estimates, where R (see Eq. (7-4) of 
Reference 8-1) instead of R (see Eq. (8-3)) is plotted to emphasize 
the relationship between resolution and w c . Clearly, the choice of T 
according to Guideline 2 provides beneficial results. 


Because Guideline 2 was developed based essentially on analysis of 
the two-sinusoid case, the second example considers its application to a 


case involving other than two sinusoids. Letting N^ = 128, N = 30, 


cos(0.7irt) + 2 cos(0.9irt) + 1.5 cos(irt) 











T is chosen so that a> c = (0.85n)T varies from O.Itt to 0.8ir. The 
results of spectrum analysis are presented in Figures 8-2a through 8-2h, 
demonstrating the usefulness of Guideline 2 when n is an even integer 
other than four. 

A common practice used to improve spectrum estimation is the 
assumption that n is larger than its actual value. The reason behind 
this is that the additional parameters estimated tend to model the error 
in the autocorrelation estimates. Although n was assumed equal to its 
actual value in the analysis leading to Guideline 2, that Guideline may 
be applied with success when n is assumed greater than its actual value. 
The third example demonstrates this. With N<j = 256, N = 100, 
a w 2 =0.1, and 

x c (t) = 2 cos(0.7irt) + cos(0.8nt) + 1.5 cos(0.9irt) 

+ 2 cos(0.95nt) + 1.5 cos(nt) 

spectrum estimates are calculated for values of T such that 
m c = (0.85ir)T varies from O.ln to 0.825ir. It is assumed that n = 20 
when spectrum analysis is performed (although actually n = 10). The 
estimates obtained are presented in Figures 8-3a through 8-3i. Although 
the results obtained when w c = ir/2 might not be considered the best, 
each spectral magnitude peak is clearly resolved, which is not the 
situation for all the values of u c considered. 

Thus, it appears that spectral resolution may be improved 
considerably if the autocorrelation sampling period is selected to place 
the autocorrelation sequence poles in desirable positions from a 
parameter sensitivity viewpoint. In addition, the results presented her 
concerning sampling period selection have added importance because of a 
result presented by Kay [Ref. 8-2J. He suggests that when calculating 
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autocorrelation estimates for a Gaussian process, the data should be 
sampled at least as often as twice the Nyqulst rate of the underlying 
process. Denoting the corresponding sampling period by T 2 N, autocorre¬ 
lation estimates would then be available at integer multiples of T 2 fj. 
Recalling that according to Guideline 2 the autocorrelation sampling 
period T should be chosen so that w c - ir/2 (assuming - 0, 
i * 1, •••,n), or 


we see that T = T 2 r when ftj is very close to fin/2 • Since the set 

{ft,, •••, ft } has been ordered such that 0 < ft, < ft,,,, i ■ 1, •••, n/2, 
1 1 n J 1 i+i 

we note that when the elements of {ftj, •••, ftn/ 2 } are tightly 
clustered (that is, when high resolution is required), setting T 
according to Guideline 2 will force the data to be sampled at a rate 
which produces good autocorrelation estimates. Thus, choosing T 
according to Guideline 2 will yield good autocorrelation estimates, in 
addition to reducing the sensitivity of spectrum analysis to errors in 
those estimates. 


The Use of Desampling to Improve Spectral Resolution 


The results presented so far concerning proper choices for T are 


obviously of no help when T is fixed or when a purely discrete-time 
process is of interest. Even if T is adjustable, it may be necessary to 
separate the poles further than afforded by constraint Eq. (8-7) to 
reduce spectrum estimation sensitivity to a reasonable level (see 
Guideline 1, p. 211 in Reference 8-1). In these circumstances it is 


necessary to desample the autocorrelation sequence to move the poles to 


desirable locations. 


In Reference 8-1, the subject of desampling was considered in 
detail, and a modified spectrum estimation procedure based on desampling 
was constructed. The essence of this procedure is the decomposition of 
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the autocorrelation sequence {r(k)| into the M desampled sequences 
{( r M,j(k)} : J = 0, • • •, M - 1} defined by 


r M>J (k) * r(M|k| + J) 

The integer M is known as the "desampling interval". Assuming the poles 
of {r(k)} are {zj, •••, z n } as in Eq. (8-1), the poles of each 
desampled sequence {r M ^(k)} will be {z^, • •*, z^}. Thus desampling 
(r(k)} has the effect of relocating each autocorrelation pole in the 
complex plane: 


If jZjJ * 1, i = l,»*»,n and the desampling interval M is chosen to place 
the poles of the desampled autocorrelation sequences in desirable 
locations according to Guidelines 1 and 2 of Reference 8-1, then we might 
suspect that spectrum analysis of the desampled sequences 
{{r M j(k)}: J ■ 0, •••, M—1} would provide higher quality results than 
spectrum analysis of the sequence (r(k)} as a whole. The modified 
spectrum estimation procedure in Section 7.4 of Reference 8-1 exploits 
this idea by estimating the parameters {A^, •••, A^} and {z^, •••, z^} 
of (r(k)} by first estimating the corresponding parameters of the 
desampled sequences. 

The improved resolution obtainable by the modified procedure when 
the desampling Interval M is chosen based on consideration of Guidelines 
1 and 2 is demonstrated here via three examples. In all three examples 
the power spectrum of a discrete-time process (x(k)} is estimated, with 
the required autocorrelation estimates being obtained from the data 
sequence {y(k) - x(k) + w(k), k » 0, •••, - l} where {w(k)} is a 

Gaussian, zero mean, white noise sequence with variance o w 2 . The auto¬ 
correlation estimator used ic defined by Eq. (8-8). Note that in light 









of Eq. (8-9), the first equation in Eq. (7-33) of Reference 8-1 as well 
as the first and (N + l)st equations of Eq. (7-18) of Reference 8-1 are 
deleted prior to obtaining fa M •••, n } (the autoregressive 
coefficient estimates for the desampled sequences) and {A^, •••, A^} for 
the examples which follow. 

The first example demonstrates the benefit of choosing M according 
to Guideline 2. We have N<j - 256, N ■ 100, a w 2 ■ 0.1, and 

x(k) - 2 cos(0.09474irk) + 2 cos(0.10530nk) 

Figures 8-4a through 8-4i contain the results of spectrum estimation for 
different values of M. Each figure contains ten spectrum estimates with 
each corresponding to a different realization of the noise process 
(w(k)}. A diagram is Included with each figure to Indicate the locations 
of the poles of |rjf,j(k)}. Noting that Guideline 2 implies that M 
should be chosen so that the poles of frjj j(k)} are "centered" around 
the imaginary axis, we see that choosing M in this manner yields good 
results. 

The second example considers the case where the poles of (r(k)) 
are already centered around the imaginary axis, but nevertheless are too 
close to be resolved (see Guideline 1). For situations such as this, 
desampling may be used to increase the angular separation of the poles 
while keeping them centered around the imaginary axis. For this example 
we have N<j » 512, N ■ 400, a w 2 - 1, and 

x(k) - 2 cos(0.49irk) + 2 cos(0.5irk) + cos(0.51irk) 

Figures 8-5a through 8-5d contain the results of spectrum estimation when 
M - 1, 13, 23, and 33, respectively. It can easily be seen that 
Increasing the angular separation of the closely spaced poles via 
desampling improves resolution considerably. 







Up to this point, all analysis and examples have dealt with cases 
where the poles of (r(k)} are on the unit circle. The final example 
demonstrates the applicability of the results of this report and of 
Reference 8-1 to those cases where the poles are slightly inside the uni 
circle. In particular, the process {x(k)> is generated by driving a 
filter of impulse response 

h(k) = 0.99^ cos ( 0>2 3Trk) u(k) 

+ 2(0.99 k ) cos(0.27irk) u(k) , — < k < +« 

with a Gaussian, zero mean, white noise process (v(k)} of variance 
a v 2 = 1 ({u(k)} is the unit step sequence). The power spectrum of 
(x(k)} is plotted in Figure 8-6. Spectrum estimates are calculated for 
various values of M, and for Nj = 512, N * 200, and a w 2 = 1. Figures 
8-7a through 8-7d contain the results of spectrum estimation, with the 
ten estimates in each figure corresponding to different realizations of 
the noise processes (w(k)} and (v(k)}. As implied by Guideline 2, 
resolution improves when M is chosen so that the poles of {r^ j(k)} are 
centered around the imaginary axis (see Figures 8-7a and 8-7b). When M 
is increased further to increase the separation of the closely spaced 
poles, small improvements in the spectrum estimates are evident (see 
Figures 8-7c and 8-7d) as implied by Guideline 1. 

8.4 Conclusion 

The results of this section together with Reference 8-1 represent 
a new approach to the problem of improving spectral resolution. Spectfi 
cally, a procedure for obtaining spectrum estimates from autocorrelation 
sequence estimates is analyzed to determine its sensitivity to errors in 
the autocorrelation estimates. The results of this analysis clearly 
demonstrate that the locations of the autocorrelation poles in the 
complex plane have a significant effect on the error sensitivity of the 
















































Power spectrum estimates of a process obtained by driving 
a filter of impulse response h(k) = 0.99^ cos (0.23 irk) u(k) 
+ (2) (0.99^) cos (0.27 trk) u(k) with unity variance, zero- 
mean, Gaussian white noise. Also, Nj = 512, N = 200, and 

o^ = l. 

w 
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spectrum analysis procedure. Specifically, it appears that error 
sensitivity becomes a particularly acute problem when poles are closely 
spaced, thus indicating the difficulty encountered when high resolution 
spectrum estimates are sought. 

As a remedy for this problem, a modified spectrum estimation 
procedure is proposed which employs autocorrelation desampling as a means 
of artificially separating closely spaced poles, thus relieving the error 
sensitivity problem. The unique features of this procedure are that it 
employs all available autocorrelation estimates, and that the size of the 
desampling interval used is not bounded by spectral aliasing 
considerations (see Eq. (7-36) of Reference 8-1). The improved 
resolution obtained by this modified procedure is demonstrated via 
numerical examples. 

Also addressed in this report is the selection of autocorrelation 
sampling periods when discrete-time techniques are used to estimate the 
spectra of continuous-time processes. It is demonstrated that if the 
sampling period is chosen to "center” the poles of the sequence obtained 
by sampling the continuous-time autocorrelation function around the 
imaginary axis, then spectral resolution can be improved substantially. 

Thus, by posing the spectral resolution problem as one of 
decreasing error sensitivity, we have observed that resolution can be 
considerably improved by placing the poles of the autocorrelation 
sequence in desirable locations either by desampling or by making an 
appropriate choice for the autocorrelation sampling period. Since 
appropriate choices for desampling intervals and sampling periods may be 
obtained easily with knowledge only of approximate pole locations, these 
results should prove to be quite valuable in situations requiring high 
resolution power spectrum estimates. 
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SECTION 9 


APPROXIMATION TECHNIQUES FOR PARAMETER 
ESTIMATION AND FEEDBACK CONTROL WITH 
DISTRIBUTED MODELS OF LARGE FLEXIBLE STRUCTURES 


9.1 Introduction 

In this section we discuss some approximation techniques that may 
be used in algorithms for parameter estimation and/or feedback control in. 
distributed models such as those arising in models typical of large 
flexible space structures. The focus of our recent efforts has been the 
development and analysis of computational algorithms, e.g., convergence 
analysis, numerical ifaplementation (software development) and testing. 
While the ideas involved are also applicable to the computation of feed¬ 
back controls, we restrict our discussions here to some of our efforts on 
techniques in the context of parameter estimation or "inverse" problems: 
given observations of a system, determine parameters in models which best 
describe structural/material properties manifested by the system in 
response to perturbations (e.g., loading). 

The importance of such problems is twofold: (i) parameter 
estimation can be viewed as a primary tool in on-orbit model development 
and analysis where one seeks to understand elastic/viscoelastic material 
properties (e.g., damping or stiffness) and to detect changes in these 
properties (e.g., due to aging or prolonged stress); (ii) parameter 
estimation is a precursor to and an integral part of development of 
sophisticated feedback control laws (via feedback operators satisfying 
infinite-dimensional Riccati equations involving functional parameters of 
the system). 
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Many of the structures of interest to aerospace engineers entail 
systems composed of composite materials in rather complex geometric/ 
structural configurations. The need for methods to investigate such 
variable-structure distributed models has, in our opinion, been clearly 
established in a number of recent efforts including References 9-1 
through 9-4. A number of investigations of parameter estimation in 
models for elastic beams has Involved approximation results (the 
Trotter-Kato theorem) from linear semigroup theory. In particular, 
problems for simple beams have been treated in this manner in References 
9-5 through 9-9. In Reference 9-10 the Trotter-Kato ideas are employed 
to establish results for hybrid models similar to those introduced later 
in this section and which are important in the study of shuttle-deployed 
payloads. However, in some instances it is advantageous to use an 
alternate approach involving a variational (weak) formulation of the 
system equations along with estimates in the spirit of those found in 
numerous papers on finite-element techniques in structural problems. In 
Reference 9-11 such a treatment was given for damped cantilevered 
Euler-Bernoulli beams. In this section we outline this approach in the 
context of models for the transverse vibration of flexible beams with 
attached tip bodies and base acceleration. 

Fundamental to our discussions is a conceptual framework in which 
one has a dynamical model of the physical system with "states" u(t,x), 

0 < t < T, x e ft, and "parameters” q(t,x), q e Q, where Q is an 
admissible class of parameter functions. The equations of motion for the 
type of problem considered generally take the form of a hybrid system of 
coupled partial (governing the vibration of the beam) and ordinary 
(describing the motion of the tip bodies) differential equations with 
appropriate geometric boundary conditions and initial data. The identi¬ 
fication problem is formulated as follows. One is given observations 
(data) uij for u(ti,Xj) and seeks to solve the optimization problem 
of finding—in the feasible parameter set Q—parameters q which give a 
best (in the least-squares sense) fit of the model to the data. By 
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feasible we mean of course that throughout the optimization, the states 
and "parameters" must jointly satisfy the dynamic constraints (i.e., the 


hybrid system of equations describing the motion of the beam and tip 


bodies). The resulting identification problems therefore, tend to have 
infinite-dimensional constraints. Moreover, if the parameters to be 
identified are functional (e.g., spatially varying stiffness and/or 
linear mass density; temporally and/or spatially varying loads), the 
admissible parameter space is a function space and as such is infinite¬ 
dimensional as well. The solution of the resulting constrained 


optimization problem, therefore, necessitates the use of some form of 


finite-dimensional approximation. We formulate this problem in an 
abstract setting with Hilbert state space V and parameter space Q. For 
computational purposes we then approximate V and 0 by finite-dimensional 
spaces V N and respectively. We illustrate these ideas with a 
specific model and particular classes of approximations in our subsequent 


The scheme we develop here is based upon the reformulation of the 
equations of motion in weak form. A cubic spline based Rayleigh-Ritz- 
Galerkin method is used to define a finite-dimensional approximation to 
the state equations. Using finite-dimensional subspaces to discretize 
the admissible parameter space, we obtain a doubly-indexed sequence of 
approximating finite-dimensional identification problems. Using standard 


variational arguments, we derive a convergence result for the state 


approximations. We show next that each of the approximating identifica¬ 
tion problems admits a solution and that from the resulting sequence of 


optimal parameter values a convergent subsequence can be extracted, the 


limit of which is a solution to the original Infinite-dimensional identi¬ 
fication problem. The approximating identification problems can be 
solved using standard techniques and readily available software. 


Numerical results which demonstrate the feasibility of our method are 
provided for a variety of examples. 











The approach described in this section represents a significant 
improvement over the method developed in Reference 9-10. Indeed, we have 
developed a scheme which is computationally simpler and, by relaxing the 
necessary hypotheses on the admissible parameter space (i.e., with regard 
to the smoothness requirements on the spatially-varying stiffness and 
linear mass density of the beam), is applicable to a wider class of 
problems. Our results are similar in spirit to those presented in 
Reference 9-12 in the context of parabolic systems, in Reference 9-13 for 
hyperbolic systems, and in Reference 9-11 for beam equations with 
standard boundary conditions (e.g., clamped, simply supported, 
cantilevered). 

We simplify our presentation by considering a cantilevered beam 
with an attached tip (point) mass. As is discussed in Section 9.4, 
however, our general approach is applicable to a broad class of beam/tip- 
body vibration problems. In Section 9.2 we derive the weak form of the 
equations of motion, define weak and strong solutions and formulate the 
identification problem. In Section 9.3 we define the approximation 
scheme and discuss convergence. Numerical results are presented in 
Section 9.4. 

We employ standard notation throughout. The Sobolev spaces of 
real-valued functions on the closed interval [a,b] whose k-th derivatives 
are L 2 (i.e., Lebesque-measurable and square-integrable) are denoted by 
H k (a,b), k - 0,1,2,... The corresponding Sobolev inner products and 
their induced norms on the spaces H k (a,b) are denoted by and 

II *11^, respectively. For Z a normed linear space with norm D • II z and 
given a map f: [0,T] + Z, we say that f e L 2 ([0,T], Z) if 

/ T ||f(t)l 2 dt < +®. Similarly, f is said to be an element in 
0 

C^([0,T], Z) if the map t f(t) from [0,T] into Z is l times 
continuously differentiable on the open interval (0,T). Finally, for a 
function of one or more real variables, the symbol D Q f (D^f) will be used 
to denote the 1st (k-th) derivative of f with respect to the independent 
variable 9. If f is a function of a single variable only, the subscript 
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may be deleted. On occasion the short-hand notation DQf(6g) or 

Df (0_) will be used In place of D.f L or DfL to denote the derivative 
0 0 |0 Q |9 0 

of f evaluated at 9g. 

9.2 The Identification Problem 
9.2.1 A Distributed Parameter Model 

We consider a long slender beam of length l with spatially varying 
stiffness El and linear mass density p which Is clamped at one end and 
free at the other end with an attached tip mass of magnitude m (see 
Figure 9-1). 



Figure 9-1. Clamped beam with tip mass; deformed. 


Using the Euler-Bernoulll theory to describe the transverse vibrations of 
the beam, we obtain the partial differential equation [Refs. 9-14, 9-15J 

p(x)D^u(t,x) + D^(EI(x) D^u(t ,x) ) «= D x (o(t,x) D x u(t,x)) + f(t,x), 

(9-1) 

x e (0,1), t e (0,T) 

for the transverse displacement u, where a denotes the internal tension 
(as a result of axial loading), and f Is the net externally-applied 
transverse or lateral load. Use of the principles of elementary 










Newtonian mechanics (i.e., force and moment balance equations) yields the 
boundary conditions at the free end. From translational equilibrium in 
the in-plane direction transverse to the beam* we have 


mD^u(t,fc) - D x (EI() 0 D^u(t, £) ) = -o(t,£)D x u(t, l) + g(t) , t e (0,T) 

(9-2) 


where g is the net external force (directed transversely) on the tip 
mass. In a similar manner, requiring rotational equilibrium about an 
axis normal to the plane of Fig. 9-1, we have 

D x u(t,*) =0 , t e (0,T) (9-3) 

The geometric boundary conditions (zero displacement and zero slope) at 
the clamped end are given by 

u(t,0) - 0 , t e (0,T) (9-4a) 

and 

D x u(t,0) =0 , t e (0,T) (9-4b) 

respectively. The initial conditions are in the form of initial 
displacement 


u(0,x) = <J>(x), x e [0, £] (9-5a) 

and initial velocity 

D t u(0,x) - «|)(x), x e [0,£j (9-5b) 

In order to characterize solutions to the hybrid system Eqs. (9-1) 
through (9-5) of ordinary and partial differential equations together 







with boundary and initial conditions, we formally represent it as an 
abstract second-order system. We make the standing assumptions that 
m > 0, EI,p e LJ0,£) with EI.p > 0, a e L 2 ([0,T], H^O,*)), 


2^ 


iKO specified in R (the set of all real numbers). Define the Hilbert 
0/, 


space H = R x h (0,1) with inner product 

< (n, d>) , (£,<[>) > R * nc + < 4>,> n 


We then rewrite Eqs. (9-1) through (9-5) as 


M Q D t u(t) + A Q u(t) = B Q (t)u(t) + F Q (t) , t e (0,T) 


Y u(t) “ 0 , Y,u(t) - 0 , Y,u(t)| = 0 , t e [0,T] 

x=0 x=0 1 |x-* (9 _ 7) 


D^u(0) 


(9-8) 


where u(t) = (u(t,£), u(t,»)) e H, F Q (t) = (g(t), f(t,«))» $ = (♦(£), <fi), 
$ = (i|>()0, i|>) and the operators M Q , A Q , B Q (t) and i « 0,1,2 on H are 


defined formally by 


M Q (n,4>) ■ (mn, P4>) 


AqU,*) - (- d(eiU) d z *(«), d 2 (ei d 2 <j>)) 


B q ( t)(Ti,4>) - (- o(t,£)D<KO, D (a D<») 

Y^TI,*) - D 1 * , i - 0,1,2 


respectively. 







There exist several ways in which the notion of a solution to the 
system Eqs. (9-6) through (9-8) can be made precise. Of particular 
Interest to us here are the ideas of a weak or variational solution and a 
strong solution. 

9.2.2 Weak Solutions 

Define the Hilbert space {v, by 

V - {(n,<|>) e H: * e H 2 (0,t), <K0) = D<|>(0) - 0, n = <f>( l )} 

< (♦(*),♦), (♦(»),*) > v - < D 2 ^. D 2 ^ > Q 

It is not difficult to show that V can be densely embedded in H. 

Choosing H as our pivot space, we have therefore the continuous 
embeddings VC H C V' , where V' denotes the space of continuous linear 
functionals on V. We now rewrite Eqs. (9-6) through (9-8) in variational 
form as the second-order initial value problem 

< M 0 D 2 u(t), 0 > H + a(u(t), 0) = b(t)(u(t), 8) + < F Q (t), 0 > H , 

(9-9) 

0 e V, t e (0,T) 

u(0) = j , D t u(0) - ijl (9-10) 


where 0 = (0(£),0) and the bilinear forms 
are giver, by 


VxV + R and b(t): VxV ->■ R 








b(t)($,$) = < -a D*, Dip > 


respectively. A solution u to Eqs. (9-9) and (9-10) with u(t) e V is 
known as a weak or variational solution to Eqs. (9-6) through (9-8). 
Indeed, if the derivatives in Eqs. (9-6) and (9-7) are taken in the 
distributional sense [Ref. 9-16], Aq and Bg(t) become bounded linear 
operators from V into V' with 

< A Q $, p > H = a($,ij) 

and 

< B Q (t)$, ip > H = b(t)($,p) 


where the H-inner-product is interpreted as its natural extension to the 
duality pairing between V and V* [Refs. 9-17, 9-18, 9-19], Since 
Fq(T) e H C V' we have therefore that the systems Eqs. (9-6), (9-7), 

(9-8) and Eqs. (9-9), (9-10) are two representations for the same initial 
value problem in V'. 

Under the standing assumptions which we have made above, standard 
arguments [Refs. 9-20, 9-21] can be used to demonstrate the existence of 
a unique solution u to Eqs. (9-9), (9-10) with u e C([0,T], h), 

D t u e C([0,T], H) and D^u e L 2 ([0,T], V’). 


9.2.3 Strong Solutions 

In order to characterize strong solutions we rewrite Eqs. (9-6) 
through (9-8) as an equivalent abstract first order system and then rely 
upon the theory of semigroups and evolution operators. Define Z = V * H 
with Inner product defined by 








We assume that El e H 2 (0,£) and o e (^([O.TJ, H^O,*)) and define 
the operator A: Dora(A)CLZ + Z by 


Dom(A) » Doib(Aq) x V 


A 


A 



I 


0 


where I is the identity on V, Mq and Aq are as they were defined 
above, and 


Dom(A Q ) - {$ = (<J>U), <j>) e V: 4 e H 4 (0,i), D 2 «(>( &) = 0} 


Similarly, define the operators B(t): Z + Z by 


B(t) 


"o‘»0 (t) 


and define A(t): Dom(A)C Z ♦ Z by A(t) » A + B(t). Let F(t) e Z be 
defined by F(t) = (0, M^^F^(t)), z^ e Z by z^ = ($,$) and consider the 
initial value problem 

D fc z(t) = A(t)z(t) + F(t) , t e (0,T) (9-11] 


z(0) = z Q (9-12] 

It is not difficult to argue that the operator A is densely defined and 
conservative; that is: 
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< Az, z >. 


z e Dom(A) 


(9-13) 


a 


Moreover, we have the following result. 

Theorem 9-1. The operator A: Dom(A)CTZ + Z is skew self adjoint. 


Proof . We first argue that -ACA*; that is: Dom(A) CZ Dom(A*) and for 
each z e Dorn (A) , A*z = -Az. Let e Do® (A) with z^ Then 


< El D iji. , D <(>„ > n - < E! D ♦. , D ♦, > 0 


~ < $1 » “ A q^ 2 >H " a($ l’ $ 2 ) 


- < z , Az >, 


1 » 2 Z 


where we have used integration by parts, the definition of Dom(AQ) and 
the definition of V in performing the above computation. This, of 
course, implies that Z 2 e Dom(A*) and A*Z 2 = -Az 2 » 

We next argue that Dom(A*)C Dom(A). Let w e Dom(A*) and set 
y = A*w. Then for z e Dom(A): 


< z,y > = < z, A*w > 7 = < Az, w >_ 


Recalling that z,w,y e Z, we can write z h (zj,z 2 ), w = (Wj,^) and 
y = (yj,^). where y ? = (yl,y|) e H and *2 “ ^2 e H * 111611 


i 

's' 


• v > 







- + < M 0 2 2> y 2 > R - aC^.Wj) + < A^, 0 2 > R 

■ < El D 2 zi, D 2 yi >o + mz2<fc)y2 + < P z 2» >0 “ 

< El D 2 z 2 , D 2 Wj > q - D(EI(£) D 2 z 1 (£))w^ + < D 2 (EI D^), w 2 > Q 

( 9 - 

Let 0 e H 2 (0,O be defined by 

D 2 0 - py 2 ; 0(0 - 0 , D0( O - - my| 

Then substituting into Eq. (9-14) and integrating by parts, we obtain 

o - - d(ei( 0 0 2 z l (O)(y 1 (O +wj) 

+ < D 2 (EI D 2 Zi ), (y x + w 2 ) > 0 + < D 2 z 2 , (0 - El D 2 Wj) > ( 

which implies that: 

(i) -D(EIU) D 2 z 1 (£))(y 1 (£) + w*) + < D 2 (EI D 2 Zi >, ( 7j + 0 2 ) > Q - 

and 

(ii) < D 2 z 2 , (0 - El D 2 Wl ) > 0 - 0 

Let (n,$) be an arbitrary element in H, Choose z\ e H*(0,O which 




satisfies 










with boundary conditions 


d(ei(0 d 2 Zi U)) = -n , d 2 Zi U) = 0 

D Zl (0) = 0 , Z j(0) = 0 

Then z^ = ( Z ^(£), ze DohKAq) and condition (1) therefore implies 
that 

< (CyjCfi.) + w*), (y L + w 2 )), (n,<l>) > H = 0 , (n,<J>) e h 

*1 „2 0 

from which it follows that yj(i) = -Wj ani * y^ “ -w 2 e H (0,0. We 

have therefore, that = -y^ e V. 

Next, let <j> be an arbitrary element in H°(0,£). Choosing 

z, - J J 4> , condition (ii) implies that 
1 0 0 


< (9 - El D 2 Wj), <j) > Q =* 0 , * e H°(0,O 


2 4 

and hence that 9 = El D w^. This in turn implies that Wj e H (0,0 and 

2 1 « 

D Wj(£) = ^^ 0(0 = 0. Since w e Z, we have w^ e V and therefore 

Wj e Dom(Ag); we conclude that w = (w^w^) e Dom(Ag) x V = Dom(A) and 

the theorem is proved. □ 


Since A is densely defined, it follows [Ref. 9-22, Theorem 3, p. 142] 
that A* is closed and, by Theorem 9-1 above, that A is closed as well. 
This fact together with Eq. (9-13) and Theorem 9-1 yield that A is 
maximal dissipative and hence [Ref. 9-23, Theorem 4.2, page 84] that it 
is the infinitesimal generator of a Cp-semigroup of contractions 
{S(t): t 0} on Z. It is in fact the case [see Refs. 9-23, 9-24] that 
S(t) is defined for t < 0 and that |S(t): - ® < t < +«} forms a 
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Cg-group of unitary operators on Z. Since the operators B(t), 0 < t £ T 
are bounded, it follows [Ref. 9-25, Theorem 2.3, page 132] that 
{A(t)}£ e [o,T] is a stable family of infinitesimal generators. Since 
Don(A) is independent of t, it follows [Ref. 9-25, Theorem 4.8, page 145] 
that the family {A( t)} t g generates an evolution system 

{U(t,s): 0 < s < t < T} on Z. 

Define 


in Z by 


z(t) = (u(t), v(t)) 

z(t) = U(t,0)Zg + / U(t ,s)F(s)ds 


(9-15) 


(9-16) 


The continuous function z given by Eq. (9-16) above is the unique mild 
solution to the initial value problem Eqs. (9-11), (9-12). If in addi¬ 
tion zq e Dom(A) (that is $ e Dom(Ao) and if> e V), then z is a strong 
solution to Eqs. (9-11), (9-12). Indeed z is differentiable almost 
everywhere (a.e.) on [0,T] with D t z e L 2 ((0,T], Z), satisfies Eq. (9-12), 
satisfies Eq. (9-11) for almost every t e [0,T], and is such that 
z(t) e Dom(A) a.e. on [0,T]. 

We shall call u given by Eqs. (9-15), (9-16) a strong solution to 
the initial value problem Eqs. (9-6) through (9-8). The following result 
is easily obtained. 


Theorem 9-2 . Suppose El e H^(0,i), a e C*([0,T], H*(0,O), $ e Dom(Ag) 
and iji e V. Then u given by Eqs. (9-15), (9-16) is the unique strong 
solution to the initial value problem Eqs. (9-6) through (9-8). We have 
that u satisfies Eq. (9-8) and Eqs. (9-6), (9-7) a.e. on [0,T]. More¬ 
over, u is twice differentiable in H and differentiable in V a.e. 
on [0,T] with D^G e L 2 ([0,T], H) and D fc G e L 2 ([0,T], V). 
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It is also not difficult to show that if a strong solution u to 
Eqs. (9-6) through (9-8) exists, then it coincides with the weak 
solution, and in such a case, it is given either by Eqs. (9-15), (9-16) 
or as the solution to the initial value problem Eqs. (9-9), (9-10). 

9.2.4 Formulation of the Identification Problem 

For ease of exposition in formulating the identification problem, 
we consider a (reasonably broad) class of inverse problems which are of 
particular interest in the structural modeling of large flexible 
spacecraft and shuttle-attached payloads [see Refs. 9-15, 9-26]. We 
assume that we are only interested in estimating the parameters m e R, 
EI,p e 1^(0,O and a e 1*2C1 0 *T] , H^O,*)) where a - a(ag,m,p) 
with ag e 1,2(0,T). It is not difficult to further generalize the 
results which follow to allow for the identification of initial data, the 
external loads f and g, and more general forms of the internal tension a 
[see Ref. 9-12]. The motivation for choosing o to be dependent upon a 
time-varying function ag, the magnitude m of the tip mass, and the 
linear mass density p of the beam is be made clear below. 

Form the Cartesian product 2= R x 1,^(0, £) x 1,^(0, i) x 1^(0,T) 
with the usual product topology. Let Q be a subset of £ which 
satisfies: 

(HI) 0 is compact; 

(H2) There exist constants , M^, i * 1,2,3 such that 

0 < ^ < m < Mj 
0 < m 2 < El < M 2 
0 < m 3 < p < M 3 
for all q - (m,EI,p,ag) e 0; and 

(H3) For all q e 0, o(q) e L2([0,T]), ^(0,0) with the mapping 
q + a(q) continuous from 0 into L 2 ([ 0 ,T], H^Ojf)). 

We assume that we ha^e been provided with displacement measurements 
{u(t.,x,): !■!,•••,y, j*l,»»»,v} for the beam at positions x e [0,£], 







j - 1,2,...v, and at times tj_ e [0,TJ, i - 1,2,.. .y, which result from 
a known input disturbance applied to the system in a known Initial 
state. The identification problem is formulated as a least-squares fit 


Find q * (m,EI,p,ag) e Q which minimizes 


J(qi u(q)) -ll lu(t ,x ) - u(t ,x ; q)H 
i-1 j-1 2 2 


subject to u(t; q) = (u(t,£; q), u(t,»; q)) being the solution to 
Eqs. (9-9), (9-10) corresponding to q e Q. 

The infinite dimensionality of both the state, which is governed by 
the system Eqs. (9-9), (9-10), and of the admissible parameter space Q 
(being a function space) necessitates the use of some form of 
approximation in solving Problem (ID). We develop and analyze one 
particular scheme in the next section. We note that the approximation 
theory to be developed below will also permit the formulation of the 
identification problem based upon criteria other than displacement; for 
example, based upon velocity [see Ref. 9-10]. 

9.3 An Approximation Scheme 

9.3.1 Approximating Identification Problems 

Our approximation scheme is based upon the construction of a 
sequence of approximating Identification problems, in each of which both 
the state constraint and the admissible parameter space are finite- 
dimensional. We argue that each of the approximating problems admits a 
solution. The resulting sequence is shown to have a sub-sequential limit 









which is a solution to the original identification problem. The state 
approximation is constructed using a spline-based Galerkin method. The 
admissible parameter space is discretized using splines as well. We 
begin by discussing the state approximation. 

Working abstractly at first, for each N « 1,2,»»», let V N CV be 
a finite-dimensional subspace of H and let P denote the orthogonal 
projection of H onto V N with respect to the inner product on H. The 
Galerkin equations corresponding to Eqs. (9-9), (9-10) are given by 

< M Q D^ N (t), 0 N > H + a(u*(t), 0 N ) = b(t)(i N (t), 0 N ) + < F 0 (t), 0 N > R 

(9-18 

0 N e V N , t e (0,T) 

i N (0) - P N | D t u N (0) - P N tp (9-19 

where u N (t) - (u^(t,£), u N (t,»)) e V N . We define the following 
sequence of approximating identification problems. 

(IDN) 

Find q N = (m N ,EI N ,p N ,aN) e Q which minimizes j(q; u N (q)), where J 
is defined by Eq. (9-17) and u N (t; q) = (u N (t,i; q), u N (t,»; q)) is 
the solution to Eqs. (9-18), (9-19) corresponding to q e Q. 

Of particular interest to us here is a scheme involving the use of cubic 
spline functions. Let denote the standard cubic B-splines on 

the interval [0,t] corresponding to the uniform partition 
A N = {°, | , . £} [Ref. 9-27], Let denote the modified 

cubic B-splines which satisfy Bj(0) * DBj(O) = 0, j * 1,2, 


N + 1; 







j - 2,3 


N + 1 


Denote = (Bj(£), Bj) and define V N C V by 


The Galerkin equations Eqs. (9-18), (9-19) take the form 


mJJw N ( 0 + Ajjw N (t) = BjJ(t)w N (t) + Fjj(t) , t e (0,T) 


where for l,j = 1,2,•N + 1, the elements of the coefficient matrices 


[Mjjy = mBj(£)B^(£) + / pB^dx 

[Ajlij - MW D^HD^dx 

t B 0 (t) llj = “ Aa(t,O0Bj)(DB^dx 

[Pj(t)] t - «(t)Bj(l) + / f (t ,•)B^dx 

[| N ] t = *(1 )bJu) + / iJ)B^dx 


















[»? N ] ± = + / <|>B^dx 

[w”]^ = B^(A)Bj(£) + / bJbJcIx 


and u N (t) has the representation 


* 

u 


(t) 


l w^(t)B I 

j=l J 


9.3.2 Convergence Results 

Our convergence arguments are based upon the approximation 
properties of spline functions. Denote S^(A N ) = span{Bj}j^| and 
S^(A N ) = spanfBj. For each function $ defined on the interval [0,£] 
let I N <(> denote that element in S^(A N ) which satisfies the interpolatory 

constraints j - 0,1.N, D(l%)(li) * 

j ■ 0,N, and let I N <j> denote that element in s3(a n ) which satisfies 
the interpolatory constraints (l^)(-jjp) ■ J = 1,2,...,N, 

D(l N <(i)(£) = D<t>(l). The interpolatory spline 1% will be well defined 
whenever <}> is well defined at the node points and D<J> at the end points. 
A similar statement can be made for I N <f>. 

We require the following two standard results concerning the 
approximation properties of interpolatory splines [Ref. 9-28]. 

Proposition 9-3 : For each <|> e h2(0,£): 

HD k (I% - *)» 0 < N" 2+k HD 2 *H Q , k = 0,1 


where C, is independent of d> and N. 









Proposition 9-4 ; For each <(> e H^(0,£): 


D k (I N <^ - N _4+k |D 4 *» 0 , k «• 0,1,2 


where C, Is Independent of (j> and N. 


The following lemma is key to establishing the desired convergence 
results. 


Then <J> ♦ 4> in H (0,1) and consequently $ -*• $ in V; 

(2) P N > I strongly In H. 

Proof . We prove assertion (1) in detail. Using the Schmidt Inequality 
[Ref. 9-28) together with Proposition 9-3, we obtain the following 
estimates: 

h* N - <Mq <. b $ N ~ <, *Ci <f>) - 

- m N $ - <m q » ni N 4> - <M 0 

_< cl. N _2 BD 2 <J>ll n 0 as N f • 


where (I N <t>) “ ((I N $)(£), 1%); and 

iid( d> N - )n 0 < id(* n - T%)» 0 + - <»>)I 0 










< kjN H$ N ~ A« 0 + c{n - 1 »D 2 <|>II 0 

< kjN »<j> N - <M q + kjN II* 1 * - <M 0 + cJn _1 BD 2 <M 0 

< k^N^HD 2 *!^ + k^'W^lg + C Jn “ 1 BD 2 4> If Q 

= ^k^J + c})N _ 1 IID 2 <(.II 0 -► 0 as N -► » 

Using the Schmidt inequality together with the first integral 
relation [Ref. 9-28], we obtain: 

iidVii 2 < 2 ii d 2 (<f> N - T n )bq + 2 hd 2 T n <m 2 

< 2k 2 N 4 n<J, N - I N 4lJ + 2 IID 2 <Mq - 2 IID 2 (d> - A)! 2 

A n 2 4—N 2 22 

< 2k 2 N n<|) - *1 q + 2k 2 N III”<(> - + 2 ID <Mq 

2 2 2 

< ct IID <Mj 

where u is independent of d> and N. Choose a function 

>p e {9: 0 e H 4 (0,£), 0(0) = D9(0) = 0}. Using arguments similar to 

those used above together with Proposition 9-4, it can be shown that 

2 N - -2 4 

ND (<|i - i|i)« 0 < kN ID W Q 

where k is independent of i|i and N. It follows that 

hd 2 (4) N - 4>)n_ < nD 2 (d> N - + hd 2 U n - ii»)+ id 2 U - *)i 








< (a + 1) nD 2 (4» - 4>) B Q + kN -2 nD%l! 0 

Since H^(0,£.) is a dense subset of h 2(0,£), can be chosen so as to 
make the first term on the right arbitrarily small and therefore that 

nD 2 (4> N - 0)l o + 0 as N » 

This proves assertion (1) of the lemma. 

Assertion (2) follows from assertion (1) and the facts that V is 
dense in H and that the P N , being orthogonal projections, are uniformly 
bounded. □ 

Theorem 9-6 . Suppose that |q N }dO with q N -► q as N + ». Let u denote 
the solution to Eqs. (9-9), (9-10) corresponding to q and let u N denote 
the solution to Eqs. (9-18), (9-19) corresponding to q . Suppose further 
that u is a strong solution (see Theorem 9-2). Then u N (q N ) -*■ u(q) in V 
and D t u N (q N ) -► D fc u(q) in H as N + <*> for each t e [0,T]. 

Proof . Our argument is similar in spirit to the one used by Dupont 
[Ref. 9-29] to demonstrate the convergence of a Galerkin method for the 
integration of a class of hyperbolic systems. We adopt the convention 
that '.e superscript N on a form or operator indicates that it be 
computed with respect to q N = (m N , EI N , p^, a^), while no superscript 
indicates that it be computed with respect to q * (m,EI,p,ap). 

Applying the triangle inequality, we have 

Hu N - ull v < llu N - P^ull v + tlpNu - u« v (9-22) 

and 

IIDijN - D^uB,. < ID - D P N ull„ + IIP N D„u - D^u II „ (9-23) 

C CH — C CH t tH 










Lemma 9-5 and the fact that u has been assumed to be a strong solution 
imply that IIP^u - ull^ > 0 and HP^D^u - D t uH^ + 0 as N -*• « for each 
t e [0,T]. 

Since u N satisfies Eqs. (9-18), (9-19) and u satisfies Eqs. (9-9), 
(9-10), for each 0 N e V N we have 


+ a N (u - P N u, 0 N ) + a(u, 0 N ) - a N (u, 0 N ) 

+ b N (t)(u N - pNu, 0 N ) + b N (t)(P N u - u, 0 N ) 


+ b N (t)(u, 0 N ) - b(t)(u, 0 N ) 


Denote u = (u (£), u M ) » P N u, denote v N = (v N (i), v N ) * u N - u , and 
N N N N 

choose 0^ ■ D.v^ e V^. Then 


+ a N (u - P N u, D fc v N ) + a(u, D t v N ) - a N (u, D t v N ) 
+ b N (t)(v N , D v N ) + b N (t)(pN»I - u, D v N ) 


+ b N (t)(G, D v^) - b(t)(u, D v^) 


1 „ W N_ -N „ -N v , N,»N ■'Nv ■> 

2 "oV • V * a <* • V >) 













- < mN(I - pN) D 2^, D^vN > H + < (M 0 - Mq)D^u, D t v N > R 

+ D t a N ((I - pN)u, vN) - a N ((I - P N )D t u, v N ) 

+ D t (a(u, v^) - a^(u, v^)) - (a(D t u, v^) - a^D^u, v^)) 

+ b N (t)(v N , D t v N ) + b N (t)(pNu - u, D t v N ) 

+ b N (t)(u, D t v N ) - b(t)(u, D t v N ) 

Integrating both sides of the above expression from 0 to t, invoking 
hypotheses (HI) through (H3) on 0, and using standard estimates we obtain 

min (nij, m 3 ) (#D t v N ll R + flv N Hy) 

< ^{maxCM^, M^)n(I - P N )D g ufl R + 

+ max(Urn - mV, Up - p N H £)• IDgUllH + »D 8 v N Hh 
+ M2#(I - pN)D g G||2 + ,vN,2 + n EI _ EX N n2. # D s ul|2 
+ llv N l|2 + C Q llCT<q N ) n nv N n 2 + C^aCq 1 *)^• IID s v n H 2 
+ c2|lcr(q N )l|2.|(i - pN)u«2 + IID s v N l(2 
+ C^Hcr(q) - o(q N )|2.uun2 + »D 8 V N I|2} ds 

+ n(i - p N )uiiJ + c 2 «v N iJ + a(i - P N )$nJ 

+ ±r- DEI - EiV.iiGnJ + c 3 «v N nJ + hei - EiV.n$nJ 


where Cq and are constants which are independent of N and C 2 and 
C 3 may be chosen arbitrarily. Choosing C 2 and C 3 sufficiently 
small (i.e, C 2 + C 3 < min (mi, m 2 , m 3 )), we obtain 



ID t v N |2 + »v N l|2 < A N (t) + / Kg( 8 ){lD 8 v N »2 + |v N lJ} ds 

where 

A N (t) - KjfKl - P^)u(t) H 2 + B (i - P N){,2 

+ I El - EI n #£- flu(t) B 2 + |EI - EI n I|2.||$|2 

+ / {II (X - P n )d 2 u(s)I |2 + max( Im - m N « 2 , up - P N l£)•BD 2 u(s )|2 

+ 11(1 - P N )D u(s) 1 2 + nEI - EI n « 2.|D u(s)l|2 
s V 00 s V 

+ Hct(s,»; q N )l|2.|l(I - P N )u(s)l|2 
+ Ba(s, •; q) - o(s,*; q N ) • n ^llu(s) |2} ds} 

and 

Kjj(s) - K 2 {K 3 + lo(s,*; q N )} 

where K^, i = 1,2,3 are constants Independent of N. Using the facts 
that q N + q In Q as B ♦ ®, that u(q) is a strong solution, that 
Lemma 9-5 holds, and that Q satisfies hypotheses (HI) and (H3), we have 
that A N -*■ 0 and that Is uniformly bounded In L 2 ( 0 ,T) as B + », 

An application of the Gronwall inequality, therefore, yields 


■ D t vN,2 + Hv n B2 >0 as N + ® 

for each t e [0,T]. Consequently, HD t u N - D t pN uH H 0 and 

lu N - P^uII v + 0 as B + ®, Together with Eqs. (9-22), (9-23), this 
completes the proof of the theorem. □ 

The continuous dependence results given in Theorems 9-7 and 9-8 
below can be verified using arguments similar to those used In the proof 
of Theorem 9-6. 







Theorem 9-7 . Let {q^fC. & with q^ ■+• q* as k ♦ <=>>. If for each N fixed, 
u N (q) denotes the solution to Eqs. (9-18), (9-19) corresponding to q, 
then u N (q^) u N (q*) in V and D^^q^) D t u N (q*) in H for each 
t e [0,T] as k ♦ ». That is, the mapping q -*■ (u N (q), D u N (q)) from Q 


into V x H is continuous for each N = 1,2,... 


Theorem 9-8 . Let {q^JCI 31 with q^ + q* as k ® and let u(q) denote the 
solution to Eqs. (9-9), (9-10) corresponding to q. If u(q*) is a strong 
solution, then u(q^) + u(q*) in V and D t u(q^) ♦ D fc u(q*) in H for each 
t e [0,T] as k + «. The mapping q (u(q), D t u(q)) is continuous as a 
mapping from 0 into V x H in neighborhoods of those q e Q for which u(q) 
is a strong solution. 


We are now prepared to prove our first major convergence result. 

Theorem 9-9 . For each N - 1,2,... fixed, problem (IDN) has a solution 

q N . The sequence {q w } admits a convergent subsequence {q^} with 
N, 

k * T e 0 as k -*• <*>. If u(TT) is a strong solution to Eqs. (9-9), 

(9-10), then q is a rolution to Problem (ID). 

Proof . Theorem 9-7 implies that the mapping q + j(q; u N (q)) from Q 
into R is continuous. This together with the fact that Q is compact 
yields the existence of a solution q N to problem (IDN). The existence of 
a convergent subsequence also follows from the compactness of 

0. If q N k + q as k + « and u(q) is a strong solution to Eqs. (9-9), 
(9-10), then Theorem 9-6 implies that 


jfa, god) 


lira Jfa k ; u k (T k )) 









J(q; u(q)) 


£ lim J(q; u^q)) - 

k+«° 

for all q e Q. The theorem is proved. □ 

Although the state equation Eq. (9-18) in problem (IDN) is finite¬ 
dimensional, the admissible parameter space 0 is a function space and 
hence the minimization of J is over an infinite-dimensional space. We 
briefly indicate a means of overcoming this difficulty which involves the 
introduction of a second level of approximation into our scheme. 

For each M * 1,2,..., define the sets (^C ^ by ^ M (Q),where 
the mappings ./ M satisfy: 

(PI) Q ■*• 2 is continuous; and 

(P2) ^ M (q) + q as M + ■, uniformly in q, for all q e Q. 

We assume further that the sets (p have the property: 

(P3) For each q e qM, u(q) is a strong solution. 

We define the following doubly-indexed sequence of approximating 
identification problems. 

(IDNM) 

Find q M N - (m^, Eljj, pN, a Q N) e qM which minimizes J(q; u N ) subject 
to u N being the solution to Eqs. (9-18), (9-19) corresponding to q. 

Note that if the sets (p are of dimension kj^ < », then the optimiza¬ 
tion specified in Problem (IDNM) is simply over a compact subset of the 
space R^ m and is subject to finite-dimensional constraints— 
a computationally tractable problem. 

The convergence arguments now go as follows. Property (PI) of the 
mappings and the compactness of Q imply that the sets Q* 1 are com¬ 
pact as well. For each N - 1,2,«»» and each M - 1,2,»»», Problem (IDNM) 
therefore has a solution q M N . For each M - 1,2,»»« fixed, the sequence 
{q M N } N “i admits a convergent subsequence with ♦ q M 
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as k > «. Recalling property (P3) and arguing as we have in the proof 
of Theorem 9-9, we conclude that JJm * s a solution to the problem of 
minimizing J(q; u) over Q M . Since Q M ■ .^(0), there exists q M e Q 
such that - «^ M (q M ). Since {<? M f C Q and Q is compact, there exists 
a convergent subsequence {q M } of {q^} with q^ + ^ as ] ♦ «. Property 
M, - 

(P2) then implies that & J(q u ) -► Tf as j ♦ “ and hence that q„ + q as 

3 m 3 

j ® as well. Since for each q e 0 • 

J(q M ; u) < J(q; u) 

M j 

M 1 M 1 

and since 0 J = & J (Q), it follows that for each q e 0: 

j(q M ; G(?T m )) < J(^ M J(q); i(^ M J(q)) (9-24 

J J 

Taking the limit as j + « on both sides of Eq. (9-24) and recalling 
Theorem 9-8 and Property (P3), we infer that for each q e Q: 

Jfa; u(T)) < J(q? u(q)) 

from which we conclude that q is a solution to Problem (ID). In summary 
we have proved the following result: 

Theorem 9-10 : Each of the approximating identification problems 
(IDNM) has a solution q M N . From the tableau a sequence fq(k)} ca 

be extracted with q(k) + q, a solution to Problem (ID), as k ♦ ». 
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Typically, an appropriate choice for the sets are finite¬ 
dimensional spaces of interpolating linear or cubic spline functions with 
the mappings being constructed from the usual interpolation 
operators. If the elements in Q are sufficiently regular, it is not 
difficult to verify that the mappings ,0^ satisfy Properties (PI) 
through (P3) for this approximation. A detailed discussion and several 
examples of this particular choice for the sets in the context of 
Inverse problems for parabolic systems can be found in Reference 9-30. 
Similar results for identification problems involving the estimation of 
functional parameters in beam equations with simple boundary conditions 
(e.g., clamped, simply supported, cantilevered) are presented in 
References 9-8 and 9-11. 

9.4 Numerical Results 

In this subsection, we discuss numerical results for a variety of 
examples. Although the analysis presented in the preceding subsections 
is based primarily on a simple example involving a cantilevered beam with 
tip mass, only minor modifications are required so as to make our general 
approach applicable to a broad class of Inverse problems for beam 
vibration. Some of these will be considered and outlined below. 

Although one of the major features of our scheme is Its ability to 
identify or estimate functional (spatially and/or temporally varying) 
parameters, our numerical findings for this important class of problems 
are, at present, incomplete. For this reason, this section includes 
examples involving the identification of constant parameters only. Our 
results for functional parameters will appear in a forthcoming paper. 

In all of the examples below, the fits were based upon artificially- 
generated "observations”. By this we mean that so-called "true" values 
for the parameters were selected and a modal-based Galerkin method was 
used to generate the solution to the resulting system of equations from 







which the sampled displacement measurements {u(t^,Xj): i«l,•••,!», 
j«l,»»»,v} were obtained. Results for examples Involving fits based upon 
actual experimental data will be discussed elsewhere. 

The flnlte-dlmensional optimization problems were solved using an 
iterative steepest-descent routine, ZXSSQ, from the IMSL Library. The 
algorithm used is the one suggested in the papers by Levenberg 
[Ref. 9-31] and Marquardt [Ref. 9-32]. The finite-dimensional second 
order initial value problem given by Eqs. (9-20), (9-21) which has to be 
solved at each iteration to compute the value of the least-squares payoff 
functional J, approximations to its gradient, and the corresponding 
Jacobian matrix, is integrated using a variable-order Adams predictor- 
corrector method (IMSL routine DGEAR). The system of Eqs. (9-20), (9-21) 
does not in general demonstrate stiff behavior. 

The integral inner products in the generalized mass and stiffness 
matrices and the generalized Fourier coefficient vectors for the external 
loads and initial data were computed using a composite two-point Gauss- 
Legendre quadrature rule. The four-subinterval support of the cubic 
B-splines and their derivatives leads to 7-banded matrices and conse¬ 
quently contributes to the overall computational efficiency of the 
method. 

All examples were run on the IBM 3081 processor at Draper Labora¬ 
tory. 

Example 1 . In this example, we consider a cantilevered beam of length 
l * 1.0 with an attached tip mass. We seek to identify the spatially- 
invariant stiffness El and linear mass density p of the beam and the 
magnitude m of the tip mass. We assume that the system is initially at 
rest C■ 0, i|» - 0), and is subsequently excited by the distributed 
lateral load along the beam given by 

f(t, x) ■ e X sin 2irt 
and the point force applied to the tip mass 







We assume that there is no axial loading (o ■ 0). 

Displacement observations were generated by taking El = 1.0, p - 3. 
and m = 1.5 to be the true values of the indicated parameters. 
Measurements were taken at positions xj * 0.25j, j * 2,3,4, and at 
times t£ = 0.5i, i = 1,2,...,10, from a solution to the system Eqs. 

(9-1) through (9-5) generated using a Galerkin method and the first two 
natural modes of the unforced unaccelerated system [see Ref. 9-26]. The 
"start up" values for the steepest descent routine were EIq =0.7, 

P0 = 2.7 and mg = 0.7. The final converged values for El, p , 
and m N together with the residual sum of squares J and required CPU 
times for various values of N are given in Table 9-1 below. 


Table 9-1. Convergence results, Example 1. 


N 


-N 

P 

-* 

m 

j" 

CPU 

(minrsec) 

2 

1.00057 

3.04455 

1.48957 

0.12 x 10" 3 

-3 

0:09.19 

3 

1.00067 

3.01256 

1.49707 

0.11 X 10 J 

„ -3 

0:22.10 

4 

1.00027 

3.00922 

1.49721 

0.11 X 10 

-3 

0:57.60 

5 

1.00016 

2.98936 

1.50262 

0.11 X 10 

-3 

1:22.52 

6 

0.99912 

2.99720 

1.49952 

0.11 * 10 

2:52.76 


Example 2 . We consider the system described in the previous example. A 
in that case, we assume that it is initially at rest and then excited by 
the distributed load 









f(t,x) 


e sin 2irt 


along the beam. The point force acting on the tip mass is (modified): 


g(t) 


2 e 


-t 


In this case, we assume that the entire system is subjected to a 
time-varying base acceleration which is given by 


a 0 ( t ) 



° < t < 1.5 
t > 1.5 


(9-25) 


The internal tension resulting from the axial load [see Refs. 9-15, 9-26J 
is given by 

a(t,x) - -a Q (t)(p( l - x) + m) 


We are interested once again in estimating the stiffness El, the linear 
mass density p, and the magnitude m of the tip mass. The true values of 
the parameters were taken to be El - 1.0, p - 3.0, and m « 1.5, with the 
reference solution being generated using the first two natural modes of 
vibration for the unforced, unaccelerated system. Displacement measure¬ 
ments were taken at positions xj * 0.625 + 0.125j, j - 1,2,3 and at 

times tj "0.51, i - 1,2.10. The start up values for the iterative 

search routine were taken as Elg - 0.7, pg - 2.7, and mg ■ 1.7. 

Our results are summarized in Table 9-2. 
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Table 9-2. Convergence results, Example 2. 


N 

ei* 

-N 

P 

—N 
m 

3" 

CPU 

(minrsec) 

2 

1.00057 

3.09966 

1.47928 

-4- 

i 

o 

o 

0:18.01 

3 

1.00121 

3.06360 

1.48727 

0.18 x 10 -4 

0:35.93 

4 

1.00092 

3.04144 

1.49207 

0.19 x 10" 4 

1:26.19 

5 

1.00057 

3.03063 

t 

1.49413 

0.19 x 10' 4 

4:15.87 

6 

1.00117 

3.03186 

1.49436 

0.35 x 10 -4 

5:21.34 


We note that although strictly speaking the convergence theory 
developed In Section 9-3 requires that ag e C*, the scheme performed 
satisfactorily with the discontinuous ag given by Eq. (9-25) above. 

Example 3 . In this example we consider a free-free beam of length 
l - 1.0 with an attached tip body at each end (see Fig. 9-2). 



Figure 9-2. Free-free beam with tip body affixed to each end. 
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The tip bodies are assumed to have known mass properties which are given 
by: 


Tip Body 0 (at x = 0) : 

mg = 0.75, eg = 0.1, Sg = ir/6, Jg = 0.6 

Tip Body 1 (at x = 1): 

mi = 1.50, ci = 0.2, 6i * ir/3, Ji = 0.4 


where for tip body i, i = 0,1, the indicated parameters represent, 
respectively, its mass m^, the distance C£ from its center of mass to 
the tip of the beam, its mass center offset 6^ as measured from the 
extension of the longitudinal axis of the beam, and its moment of inertia 
about its center of mass. 


The equations describing the transverse displacement of the beam and 
the translational and rotational equilibrium of the tip bodies are given 
by [see Refs. 9-15, 9-26]: 


pD 2 u + D 2 (EI D 2 u) = f , x e (0, 1), t e (0,T) 


(9-26) 




s6 0 »t D x U + D x( EI D x u ) “ 0 * x - 0, t e (0,T) 
..i 0 c 0 cos6 0 D^u + (J Q + m Q CQ)D^D x u - El D 2 u = 0 , x - 0, t e (0,T) 

x = 1, t e (0,T) 


ra.D^u + me.cos6. D 2 D u + 0 (El D 2 u) = 0 
It 11 1 t x x v x ' 


mjCjCosfij D 2 u + (Jj + mjC 2 )D 2 D x u + El D^u =0 , x * 1, t e (0,T) 


The initial conditions are of the form 





V'r 

iUi 

Wei 










u - ♦ , D t u - ^ , X e [0,1], t - 0 (9-27) 

Setting H = x H°(0,1) with Inner product 

< (n,4>)> (c,<|») > H = n T ? + <4>,i|<> 0 

and 

V - {$ e H: $ = (n,<(>), ip e H 2 (0,1), n - (*(0), D<K0), $(1), D<K 1 ) ) T } 

with Inner product 

< > v - < + ,+ > 2 

the weak form of the system given by Eqs. (9-26), (9-27) above becomes 
< M Q D2G(t), 0 > H + a(u(t), 8) - < F Q (t), 0 > H , t e (0,T), 6 e V (9-28) 

u(0) “ $ , D t u(0) ■ $ (9-29) 

where the Indicated operations are defined as follows: 

u(t) - ((u(t,0), D x u(t,0), u(t,l), D x u(t,l)) T , u(t,0) e V 




M n (n,<t>) 


(M n tl, P*) 










a: V x V -*■ R is given by 


a($, |) - < El D 2 (j), D 2 t(j > q 

F Q (t) = (0, f(t,.)) 

♦ = ((<K0), D(J>(0), 4,(1), D4>(1)) T , 4>) 

and 

$ - ((4>(0), D4»(0), 4,(1), D4,(1)) T , ♦) 

The appproximation scheme for the system Eqs. (9-28), (9-29) is 
constructed in essentially the same manner as it was in Section 9-3. 

In this example we seek to identify El and p. The system is assumed 
to be initially at rest and subsequently excited via the input 
disturbance 


f(t,x) * 10 e x sin 2irt 

Observations at positions xj = 0.25j, j - 0,1,...,4 and at times 
t£ » 0.21, i « 1,2,...,10, were generated using the first six natural 
modes of the system: rigid-body translation (transverse to the beam—in 
plane) and rotation (about an axis normal to the plane) and four flexible 
modes. The true values of the parameters were as; umed to be El * 1.0 and 
p - 3.0 with start up values taken as EIq * 0.7 and pq - 2.7. Our 
results are given in Table 9-3. 
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1.0 


a n (t) S 


, 0 < t < 1.5 
0 , 1.5 < t < 3.0 
1.0 , 3.0 < t < 4.0 
0 , 4.0 < t 


The equations of motion are given by [see Ref. 9-15] 


pD^u + El D^u = -a^D (p(l - x) + mjD^u + f , x e (0,1), t e (0,T) 


mc(cos6)D^u + (J + mc^D^D^u + El D^u = mc(cos6)CgD x u + mc(sin6)aQ, 

x = 1 , t e (0,T) 

(9-30 

with boundary and Initial conditions 

u-0 , Du = 0 ; x = 0, t e [0,T] (9-31 


and 


u ■ 4> , D u - t|» ; x e [0,1], t - 0 


(9-32 


respectively, where the quantities m, c, 6, and J are as they were 
defined in Example 3. Once again only minor modifications (similar to 
those outlined in the previous example) are necessary so as to make the 
theory presented in Sections 9-2 and 9-3 applicable to the system given 
by Eqs. (9-30) through (9-32) above. Taking the true values of El and p 
to be 1.0 and 2.0, respectively, and setting m * 4.0, c = 0.2, 5 * ir/3, 
and J - 0.4, displacement observations at positions x^ ■ 0.625 + 0.125J, 
j - 1,2,3 and at times tj = 3.5 + 0.51, i ■ 1,2,3 were generated using 
the first three natural modes of the unforced, unaccelerated system. 


S 

I 

£ 













The start-up values were taken as Elg * 0.7 and pg = 2.5. The final 
converged values for EI N and p N together with the residual sum of squares 
J N for several values of N are given in Table 9-4 below. 

Table 9-4. Convergence results. Example 4. 


N 

El* 

-N 

P 

j* 

2 

1.00046 

2.07731 

0.40 x 10~ 3 

3 

1.00135 

2.06652 

0.22 x 10 -3 

4 

1.00117 

2.05180 

0.23 x 10 -3 

5 

1.00086 

2.04776 

0.26 x 10~ 3 
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SECTION 10 


LARGE ANGLE SPACECRAFT MANEUVERS USING AN ANALYTIC SOLUTION FOR 
TRACKING AND DISTURBANCE ACCOMMODATION 


10.1 Overview 

In this section we consider the problem of maneuvering a 
flexible spacecraft through a large angle, where the disturbance- 
accommodating feedback control tracks a desired output state. The 
desired output state Is provided from an open-loop solution for the 
linear system model. The components of the disturbance vector are 
assumed to be represented in terms of Fourier series. Closed-form 
solutions are provided for the Riccati, prefliter, state trajectory, and 
residual-state trajectory equations which define the optimal control. 
Example maneuvers are presented where control-rate penalties have been 
Included In the performance Index for frequency-shaping, In order to 
smooth both the open-loop and closed-loop control commands. 

A common method for investigating the effectiveness of various 
control designs consists of studying frequency domain characteristics of 
the control, by numerically evaluating the required Fourier transform. 
For finite-time open-loop and closed-loop control problems this can be 
accomplished by either numerically Integrating the Integral definition o 
the Fourier transform for each frequency of Interest, or by using a 
fast-Fourler-transform algorithm. Alternatively, we present In this 
section a computationally efficient closed-form solution for the Fourier 
transform of finite-time open-loop and closed-loop control problems, 
where the dynamics of the control Is governed by matrix exponentials. 







10.2 Introduction 


The determination of a control time history for maneuvering a 
flexible spacecraft has been the subject of many recent investigations In 
the area of spacecraft attitude control [Refs. 10-1 through 10-24], Many 
early papers developed open-loop methods [Refs. 10-1 through 10-11, 

10-13, 10-19, 10-24] by assuming that the control can be obtained from 
either Pontryagin's principle or by assuming that the control can be 
expressed in terms of a Fourier series. More recent papers have 
developed closed-loop methods [Refs. 10-6, 10-10, 10-12, 10-15 through 
10-18, 10-20 through 10-23] which consider classical regulator, tracking, 
and terminal-tracking control formulations. Moreover, in an effort to 
minimize the deleterious effects of control spillover, a number of 
authors have proposed frequency-shaping techniques [Refs. 10-2, 10-14 
through 10-17, 10-19, 10-21 through 10-25, 10-27], which attempt to 
roll-off the high frequency content of the applied control. 

The dual tracking and disturbance-accommodation control problem 
described in this section extends the closed-loop tracking results of 
Turner, Chun, and Juang [Refs. 10-15, 10-16] by permitting simultaneous 
frequency-shaped tracking and deterministic disturbance rejection. The 
resulting control algorithm is applied to the problem of maneuvering a 
flexible spacecraft through large angles about a single axis. 

In Section 10.3 the dual tracking and disturbance-accommodation 
control problem is formulated by defining an inhomogeneous linear- 
quadratic optimization problem, where the plant dynamics equation is 
forced by the control and a disturbance vector. The state variable 
description for the disturbance is specified in terms of a Fourier 
series. The performance index for the control problem consists of 
integral quadratic terms in the state, control, and control-rates [Refs. 
10-16 through 10-18, 10-21 through 10-23, 10-25, 10-27]. The necessary 
conditions defining the solution for the control problem are given by two 
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coupled time-varying Riccatl-like differential equations [Ref. 10-28]. 

The first equation Is the standard Riccati equation, which possesses a 
well-known solution In terms of a steady-state plus transient term 
[Refs. 10-12, 10-14 through 10-18, 10-28 through 10-34]. The second 
equation Is a vector prefllter equation which has foreknowledge of the 
desired output state and of the deterministic disturbance which Is to be 
suppressed. A closed-form solution is presented In Section 10.4 for the 
pre-filter equation. In Section 10.4.2 the solutions for the Riccati and 
prefilter equations are cast In recursive form for efficient computer 
Implementation. Closed-form solutions are presented for the state 
trajectory equation In Section 10.5 and for the residual-state trajectory 
equation in Section 10.6. Maneuver simulations are presented In Section 
10.7 for frequency-shaped control-rate penalty designs [Refs. 10-14 
through 10-18, 10-25, 10-27]. An analytic Fourier transform for the 
control is presented in Section 10.8 for a class of finite-time control 
problems. 

10.3 The Optimal Tracking and Disturbance-Accommodation Control Problem 


The optimal reference tracking problem is formulated by finding 
the control inputs u(t) to minimize 


Y l| y*( t f ) - yO^ig + j I f {ny*(t) - y(t)«Q+ nu(t)n£ 











where IIpII w = p Wp, y is the desired output state, A is the system 
dynamics matrix, B is the control influence matrix, d is the disturbance 
state which is modeled in terms of Fourier series, C is the output 
matrix, S is the terminal state weight matrix, Q is the state weight 
matrix, R is the control weight matrix, x is the N-dimensional state 
vector, and u is the N c -dimensional control vector. Moreover, we 
assume that the system is controllable and observable. 

In order for the control system to track in the presence of known 
disturbances, we need a state variable model which generates the deter¬ 
ministic disturbance that excites the spacecraft response. To this end, 
we model the disturbance state vector as [Ref. 10-19] 

d(t) - Ae Mt 8 (10-3) 

where A is a constant matrix whose rows consist of the Fourier series 
coefficients for individual disturbances, M is a constant block-diagonal 
matrix containing the frequencies of the disturbances, and 8 is a con¬ 
stant vector. Note, however, that Eq. (10-3) is a general expression for 
either a Fourier series or a trigonometric series which models sub¬ 
harmonic terms for nonlinear effects [Ref. 10-23]. 

For maneuvers where the state is augmented by the control and 
control-rates for frequency-shaping, the A and B matrices are modified as 
shown in References 10-14 through 10-18 and 10-27. 

The desired output vector is modelled by 

y (t) - Cx (t) - CHe nt 8 0 (10-4) 

s Q - (x*(0), X*(0)) T 

where H * [I • 0] is a N x 2N selection operator, 0 is a 2N x 2N constant 
matrix, and A(t) is the open-loop costate vector [Refs. 10-8 through 
10-11, 10-13]. 
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As shown in what follows, the explicit presence of a matrix expo¬ 
nential in Eq. (10-4) permits a closed-form solution to be obtained for 
the prefilter equation of the tracking/disturbance-accommodation control 


10.4 Necessary Conditions and Solutions for the Optimal Control 
The necessary conditions defining the tracking/disturbance- 
accommodation problem for Eqs. (10-1), (10-2), and (10-3) are given by 
the following Riccati-like differential equations [Ref. 10-28]: 


P(t ) - C SC 


l - -[A - BR -1 B T P] T £ - C T QCx (t) + Pd(t) 


and the optimal control is given by 


-R~V[Px - c] 


10.4.1 Closed Form Solution for the Riccati Equation 


As shown in References 10-10, 10-14 through 10-18, and 10-27 
through 10-34, the solution for Eq. (10-5) can be written in terms of 
steady-state plus transient term, as follows: 


P(t) - P + Z (t) 
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where P 88 is the solution for the algebraic Riccati equation 
[Refs. 10-28, 10-29, 10-31, 10-33, 10-351 


T -IT T 

-PA- A P + P BR BP - C QC = 0 

ss ss ss ss 

and the closed form solution for Z(t) can be shown to be [Refs. 10-33, 
10-36, 10-37]: 

A(t-t f ) A T (t-t ) 

Z(t) ‘ Z ss + 6 l^f) - z s J e <10-1 

where A = A - BR *B T P gg is the system stability matrix, 

Z(t f ) = (C T SC - Pgg)" 1 * J * AZ + ZA T - BR -1 B T , and Z gg satisfies the 
algebraic Lyapunov equation [Refs. 10-26, 10-33, 10-38, 10-39] 

— —T — 1 T 

AZ + Z A = BR B 

ss ss 


10.4.2 Closed Form Solution for the Prefilter Equation 

The solution for £ in Eq. (10-6) follows on assuming the product 
form solution [Refs. 10-16, 10-17] 

£(t) - Z -1 (t)r(t) (10-' 

where Z is defined by Eq. (10-8) and the linear constant-coefficient 
vector differential equation for r can be shown to be 


Z(t)P d(t) + d(t) - Z(t)C QCx (t) 



r - Ar 








The solution for r follows as 


r(t) = e At [r 0 - Y (t)] 


(10-10) 


where 

y(t) = / e At {Z(x)C T QCx*(x) - d(x) - Z(x)P d(x)} dx 

0 38 

and the initial condition for r is given by 

f T * 

r Q = e Z(t f ) C A SCx (t f ) + Y (t f ) 

The solution for the integral expression in Eq. (10-10) follows on 
introducing Z(t) from Eq. (10-8), x*(t) from Eq. (10-4), and d(t) from 
Eq. (10-3) into Eq. (10-10), leading to 


Y(t) = 'f(-A, Dj, 0, t)s Q + D 2 T(A T , D 3 , 0, t)s Q 

- 4'(-A, D 4 , M, t)8 - D 2 f(A T , D 5 , M, t)B 

where 

T(A, B, C, t) = / e AT Be Cx dx 

0 

a t a "At. ~A T t f 

D = Z C 0CH , D = e [Z(t ) - Z ]e 


( 10 - 11 ) 


(10-12) 




p 

1 

§ 









and A is x B is n^ * nj, C is n^ * n£, and ¥( •) is n^ * n^. The 
matrix exponential solution for ¥(•) in Eq. (10-12) is presented in 
Section 10.9.1. Also presented there are the recursive equations for 
generating the solution for ¥(•) at discrete times. 

Subject to the recursion relationships of Section 10.9.1, the time 
histories for P and £ in Eq. (10-7) easily follow, where Z[(n + l)At] in 
Eq. (10-8) is given by 


“At, „ , -A t, 

Z(0) = Z ss + e f [Z(t f ) - Z S8 e f 


10.5 Closed-Form Solution for the State Trajectory Equation 

The closed-loop state trajectory equation is obtained by 
introducing Eq. (10-7) into Eq. (10-2), yielding 


x = (A - BR Vz" 1 ^) x + BR V5(t) + d(t) (10-13) 


The solution for x(t) in Eq. (10-13) is obtained by carrying out two 
steps. First, the homogeneous part of Eq. (10-13) is solved in terms of 
a state transition matrix solution. Second, using a variation of 
parameters technique, the complete solution for x(t) is obtained. 

As shown in References 10-10, 10-12, and 10-18, the homogeneous 
solution for Eq. (10-13) is given by 







<f>(t) = Z(t)e Z "(0) (10-14a) 

$(t) = [A - BR _1 B T Z _1 (t)]<(.(t) ; <j.(0) = I (10-14b) 

Denoting 

-T 

♦ j(t) = <f> _1 (t) = Z(0)e A Z Z _1 (t) (10-15a) 

it follows that 

*j(t) = -^(OfA - BR _1 B T z' 1 (t)] ; 4^(0) = I (10-15b) 

Equation (10-15) is recorded above since it is required for the solution 
of Eq. (10-13). 

Let us assume that the general solution for Eq. (10-13) is given 
by 

x(t) = <|>(t)[x 0 + x p (t)] (10-16) 

where x p (t) is to be determined. Introducing Eq. (10-16) into 
Eq. (10-13), while recalling Eqs. (10-14) and (10-15) leads to 

x p (t) = 4> 1 (t)BR _1 B T C(t) + ^(t)d(t) 

where the solution for £(t) is defined by Eq. (10-9) and d(t) is defined 
by Eq. (10-3). The solution for x p (t) can be shown to be 

x_(t) = d>,(t)r(t) + r,(t). + r.(t) + x__ (10-17) 















where r(t) is defined by Eq. (10-10), and 


r (t) - / 4 (t)Z(t)C T Q y*(x) dx (10-18) 

0 1 

r (t) - -/ <J> (x)Z(x)P d(x) dx (10-19) 

£ 0 

x p0 " “ r o ( See Eq * (10-10) ) 


Substituting Eq. (10-17) into Eq. (10-16) yields the following solution 
for the state trajectory equation 

x(t) = <f>(t){x Q + rj(t) + r 2 (t) + x pQ } + r(t) ( 10 - 20 ) 


The solutions for rj(t) and r 2 (t) are obtained by introducing 
Eq. (10-15) into Eqs. (10-18) and (10-19), leading to 

r x (t) - Z(0)’f(A T , D 3 , 0, t)s Q 
r 2 (t) - Z(0)'P(A T , D 5 , M, t )8 


(10-21) 


( 10 - 22 ) 


where is defined by ” 1 . (10-12). Moreover, we observe that the 

integrals appearing in Eq. .10-21) and (10-22) have been previously 
defined in Eq. (10-11); thus no additional computational effort is 
required in order to produce r^(t) and r 2 (t). 
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Assume that the residual state trajectory equation Is given by 


where u is defined by Eq. (10-7) and d r ■ A r eMtg is analogous to 
Eq. (10-3). By introducing Eq. (10-7) Into Eq. (10-23) we obtain 


where x is defined by Eq. (10-20). The solution for x can be shown to be 


A r t A r t ,t -A r t _t t. 

V‘> ‘ • «r0 - * ' / * " ” ’ ' 


f l e rT B r R 1 B T [P gg x(T) + Z~ 1 (t){x(t) - r(r)}]di 


, r A r c “ A r T Mt , 
+ {e / e A e dx}I 

0 


where r is defined by Eq. (10-10). After considerable labor, the 
integrals above can be expressed as 


x (t) = e r x _ + l 1,(0 

r ru ±ml i 


I.(t) - -[r,(t)l_a , I„(t) - -[r.(t)],,s rt , I,(t) S [r (t)],,8 














VO = -lr 3 (t)] 12 D 2a 


, I 5 (t) *-[r 3 (t)] 12 r 0 , I 6 (t) - tr 3 (t)] 13 s 0 

i 7 (t) - -[r 4 (t)i 13 6 , i 8 (t) = -tr 5 (t)l 12 a , i 9 (t) £ -[r 5 (t)] 13 s ( 

I 10 (t) " Ir 6 (t)1 13 6 » 1 ll (t) " eArt * A r» M » t)B ’ 

a 2 Z _ 1 (0)(K 0 - r 0 ) 

r^t) - r(A r , d 6 , -a t , d 3 , n, t) , r 2 (t) - r(A r , d 6 , -a t , d 5 , m, t) 

r 3 (t) - r(A r , d ? , a, Oj, n, t) , r 4 (t) = r(A r , d ? , a, d 4 , m, t) 

r 5 (t) - r(A r , D g , -a t , d 3 , r ( t) , r 6 (t) = r(A r , D g , -a t , d $ , m, t) 

and r is defined by Eq. (10-48), [r]j 2 " G^(t) is defined by Eq. (10-46) 

[r ]13 - (t) is defined by Eq.(10-47), and 

D - B R -1 B T P Z D - B R _1 B T P D = B k” 1 ® 1 

6 r ss ss 7 r ss 8 r 


10.7 Illustrative Examples 

The specific model considered in this subsection consists of a 
rigid hub with four identical elastic appendages attached symmetrically 
about the central hub, and is derived from the experimental structure of 
Reference 10-41 using NASTRAN data (see Fig. 10-1). In particular, the 
following idealizations are considered: (i) single-axis maneuvers; 

(ii) in-plane motion; (iii) anti-symmetric deformations; (iv) small 








linear flexural deformations; (v) only the linear time-invariant form of 
the equations of motion are considered; and (vi) the control actuators 
are modelled as concentrated torque-generating devices. Figure 10-2 
shows the first three antisymmetric modes which, together with the rigid 
body mode, defines the full-order model. The control system for the 
vehicle consists of a single controller in the rigid part of the 
structure. The structural parameters of the model are given in Reference 
10-41. Because of the above assumptions, only the anti-symmetric modes 
are used for the example cases of this subsection. In addition, 
full-state feedback is assumed. 

The output vector is assumed to be given by 



y^ ■ central hub angular velocity 

y 2 “ angular velocity at the mid-span location 

relative to the central hub 

y^ = central hub angular position 
y^ = tip deflection of each appendage 
y^ = controls and control-rates 

The output sub-vector y 5 is present in the formulation when control- 
rates are penalized. The elements of yj correspond to the additional 
augmented states. In the following cases, the augmented states consist 

of the control torque and control-torque rate, since the performance 

index includes a penalty on the second time derivative of the control 
torque. 
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The state vectors x and x* consist of the following elements: 



where 

x^ » central hub angle 

*2 * first mode amplitude 

x^ ■ central hub angular velocity 

x^ » first mode amplitude rate 

x^ * control torque 

x^ * control torque rate 

The second and third flexible modes are treated as unmodelled residual 
modes, which are excited by the control torque and which affect the out¬ 
put variables y \, y^, y 3 , and . 

Table 10-1 shows the weight matrices which are used for all of the 
test cases. Each test case involves a 5 second, 10° maneuver, starting 
at rest and ending at rest. The graphical results of the test case 
maneuvers are shown in Figs. 10-3 through 10-7 and are described 
qualitatively in what follows. 

Case 1 (Fig. 10-3) shows the results for a tracking controller 
which does not take into account the presence of disturbances. The 
disturbance in this case is provided at the central hub by a disturbance 
torque described by 

Tj * 0.1 sin 2irt + 0.1 sin 4irt 










diag [0.1, 

1.0 

0.5, 1.0, 0.1 

l, l(-9), l(-9)] u; 


diag [0.001 

, 1.0, 0.001, 

, 1.0, l(-9), l(-9)] 


1.0 




“5.38(9) 




1.83(9) 

1.86(9) 
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1.58(9) 
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3.63(9) 


2.22(9) 
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Figure 10-5. Case 3: Linear tracking controller with multipl* 
disturbance torques. 
























Notice that the hub angular velocity, mid-span relative angular velocity, 
and tip displacement do not reach their desired zero values. Also note 
that the first and second torque rates have spikes near the final time in 
their time histories. Undoubtedly, the spikes may be reduced by lowering 
the terminal weight matrix. However, this increases the error in the 
final states. The fact that no spikes occur in the output states and 
control-torque time histories must be attributed to the use of control- 
rate smoothing, since the large spike in the second torque rate is 
integrated twice to obtain the control torque. 

Case 2 (Fig. 10-4) presents the tracking/disturbance-accommodation 
controller which is exposed to the same disturbance as shown in Case 1. 
Since this controller utilizes disturbance information in its controller 
design, it is able to meet the desired final conditions. Note that the 
state and control time-histories are nearly symmetric about the mid- 
maneuver time, as the open-loop reference solution should be. Notice 
that the control and control-rates are much smoother than for Case 1, and 
that the peak torque requirement is slightly lower than for Case 1. 

Case 3 (Fig. 10-5) shows the simple tracking controller subjected 
to multiple disturbances, one at the central hub and one at the mid-span 
of each of the four appendages (Fig. 10-6). The appendage disturbances 
are all identical in order to satisfy the anti-symmetric deformation 
assumption. The central hub disturbance torque is given by 

T. = 0.1 + 0.1 sin 2irt + 0.2 cos 2irt 

d 

and the appendage disturbance torques are given by 

Tj = 0.01 sin 2irt + 0.02 cos 2irt 
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Notice that the output states deviate slightly from their desired final 


values. Also note the large spikes in the first and second torque 


rates. The high frequency content seen in the mid-span relative angular 


velocity time history is due to the unmodelled residual modes which are 


excited by the jump discontinuities in the disturbance torques at the 


Case 4 (Fig. 10-7) shows the tracking/disturbance-accommodation 


controller subjected to the disturbances of Fig. 10-6. In this case the 


states and control achieve their desired final values, and the control 


and control-rate time histories are smooth. Note that since the central 


hub disturbance torque has a non-zero bias and is non-symmetric, the 


resulting state and control time-histories are no longer symmetric about 
the mid-maneuver time. The residual mode response is also evident in the 


mid-span relative angular velocity time-history. Also note that the peak 
torque requirement is 10% lower than for the simple tracking controller 


We further observe that the formulation of Section 10.4.2 reduces 


to the special case of a tracking solution by setting 8 = 0 in 
Eq. (10-11), or to the special case of a disturbance-accommodation 
solution by setting sq = 0 in Eq. (10-11). 


10.8 Analytic Fourier Transform for the Control 


The following subsection provides an analytic Fourier transform 
for the control with a class of finite-time control problems. The 
fundamental definition of the complex Fourier transform is given by 


u(a>) = / dt u(t)e 


vhere u(t) is assumed to be given by [Refs. 10-9 through 10-12] 









l(t) 


. DC 
Ae c 


(10-26) 


In the above, A is n x m, B is ra x m, c is m x 1, and e^’^is the matrix 
exponential function. 

Introducing Eq. (10-26) into Eq. (10-25) yields 

u(w) - AC(a>) (10-27) 

where 

C(a>) - / T dt e Bt ce" ia)t (10-28) 

0 

The integral appearing in Eq. (10-28) can be evaluated by the 
approach given in Section 10.9.1. In particular: 



from which it follows that 


5(u>) = e BT Gj(x; m) 


(10-31) 


243 




im 










Since the numerical effort required to compute Gj(»; w) for each 
desired value of oi is prohibitive, we present in the next section a 
coordinate transformation which greatly reduces the computational burden 
required to produce £(ui) in Eq. (10-31). 
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Thus, the problem of determining $ becomes that of solving the following 
linear equation for p: 


[B - iul]p = -c (10-33) 

which leads to a well defined solution for p provided that iu> is not an 
eigenvalue of B. 

From Eq. (10-32) it follows that the matrix exponential of 
Eq. (10-30) can be written as 



Comparing Eqs. (10-30) and (10-34) it follows that the desired integral 
for £(w) in Eq. (10-31) is given by 

5(a>) = p(«j) - e BT p(a))e _1 “ T (10-35) 

where the entire solution follows after determining p(w) from 
Eq. (10-33), for each frequency of interest. The significant feature of 
Eq. (10-35) is that the computationally intensive solution for e Bx must 
be carried out only once, thus greatly reducing the labor required to 
produce ?(<*>)• 


10.8.2 Solution for the Uncoupling Transformation Vector 


Since the B matrix in Eq. (10-33) is constant and generally fully 
populated, we seek a solution technique which minimizes the computational 
effort. However, we recognize that there are two classes of solutions 
possible. First, systems where B is diagonalizable; and second, systems 
where the eigensystem for B is ill-conditioned. 
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The solution for both classes of problems are obtained via a 
“transformation method”. In particular, such methods are based upon the 
equivalence of the problems [Refs. 10-26, 10-39] 

[B - iidl]p(u») = -c 

and (10-36) 

[A - lwl]p(o)) - -y 

where A Is the appropriate canonical form of B under the appropriate 
transformation, with p(to) and y arising from the same transformation. 

The solution algorithms for both classes of problems are listed In 
Table 10-2. However, if to Is an eigenvalue of B, then Eq. (10-30) must 
be used to obtain a solution. 

The desired solution for u(w) follows upon introducing Eq. (10-35) 
into Eq. (10-27), yielding 

u(u>) = A{p(w) - e BT p(o))e~ ia>T } (10-37) 

In order to evaluate Eq. (10-37) efficiently. It is necessary to recast 
the equation in the form 

u(oj) = Ap(w) - Ap(w)e lu)T (10-38) 

where A = Ae Bt is n x m with n < m. 

10.8.3 Example Application 

Given the first-order system 

x - -x + u ; x(0) - 0 , x(t) - 1 
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Table 10-2. Solution techniques for p(u>). 


B Diagonalizable 

B Ill-Conditioned Eigensystem 

- diag (X ,...,A n ) 

A * Upper Quasi- 


triangular* 

T 

„T 

- L c 

Y - U c 

- Y /(A - iu>) , j - l,...,n 

[ A - itol] p = ~Y 

- Rp(oi) 

p(oj) - Up(w) 

« Right Eigenvector of B 

U T BU - A (real Schur 


decomposition) 

*= Left Eigenvector of B 

U T U = I (orthogonal) 

■ I (Biorthogonality) 



A quasi-triangular matrix Is triangular with possible nonzero 2x2 
blocks on the diagonal. 











we seek the control u to minimize 


J 



dt 


The open-loop optimal control can be shown to be 
u(t) - -X(t) 

where A(*) is the costate trajectory, satisfying: 



and A(0) - -l/sinh t. Thus the control of Eq. (10-39) can 

as 

u(t) * Ae Bt c ■ e t /sinh t 

where A - [0 -1], e Bt is defined by Eq. (10-40), and 

c - (0 -l/sinh t) T . The analytic Fourier transform of u(t) 


(l-iw)T 


(10-39) 


(10-40) 


be written 


follows as 


u(oj) - (/ 


(l-iw)t 


dt)/sinh 


- 1 ]/[(1 - ltd) sinh t] 
(10-41) 






Since B is diagonalizable, we use the right and left eigenvector 
transformation method to solve Eq. (10-36), leading to 



1//2 

-1//2 


21/2 

0 


R = 

0 

21/2 

, L - 

1//2 

1//2 

. A - 


i 




from which It follows that 

Y - -L T c * (l/(/2 slnh t) , l/(/T sinh t) ) T 

p (to) = C~l/[(1 + t»)^2 sinh t] , 1/ [ (1 - io))/T sinh x]) T 

p(w) - Rp(io) - 

T 

■ 1/ [2(1 + iu>) sinh t] - 1/ [ 2( 1 - iu>) sinh t] , 1/ ((1 - iu>) sinh r]) 


A = Ae BT = [ 0 -e T 

Thus, from Eq. (10-38) we have 


u(w) * [e^ ia) ^ T —1 ]/[(1 - iw) sinh t] 


We note that Eq. (10-42) agrees with Eq. (10-41). 









10.9 Appendix 


10.9.1 Closed-Form Solution for the Prefilter Equation Matrix Integrals 

In order to evaluate the prefilter equation presented in 
Section 10.4.2, we must compute four matrix integrals of the form 


'V(t) 


*F(A, B, C, t) 



(10-43) 


where the (constant) matrices A, B, and C are specified. The integral 
Involving the matrix exponentials in Eq. (10-43) can be easily evaluated 
by forming the constant matrix 


n 



»1 i" 


l 


} n 


2 


and computing the corresponding matrix exponential e^t. As shown in 
Reference 10-40, the latter has the representation: 


where 


Fj(t) Gj(t) } n { 

0 F 2 (t)J } n 2 
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Thus the solution for the integral in Eq. (10-43) follows as 


Y(A, B, C, t) = e At Gj(t) (10-44) 

In order to evaluate Eq. (10-44) at discrete time steps, the 
semi-group properties of exponential matrices are exploited, yielding the 
following recursion relationships for the matrix partitions of e nt : 

F i (t + At) = F 1 (At)F 1 (t) , F 1 (0) ■ I , i - 1, 2 

Gj(t + At) - F 1 (At)G 1 (t) + Gj(At)F 2 (t) , G^O) * 0 

where At - (t f -t Q )/m and m Is the total number of discrete time steps. 
The matrix F^(t) - e At is obtained in a separate calculation. As a 
result, the integral in Eq. (10-43) can be written (at t * (n + l)At) as 

Y(A, B, C, (n + 1)At) - F 3 [(n + 1)At] Gj(n + 1)At] 

Since each Psi-function in Eq. (10-11) is post multiplied by a vector, 
we can write 

f(A, B, C, (t + At))a = F ? (t + At)Vj(t + At) (10-45) 

where a is a vector and the recursion relationship for is given by 
v (t + At) = F (At)v (t) + G (At)v-(t) ; v (0) = 0 


I 










v 2 (t + At) « F 2 (At)v 2 (t) ; v 2 (0) = a 


A separate recursion relationship of the form of Eq. (10-45) is 
set up for each term in Eq. (10-11) consisting of the product of a 
Psi-function times a vector. 


10.9.2 Closed-Form Solutions for the Residual-State Trajectory Matrix 
Integrals 

In order to evaluate the residual-state trajectory equation pre¬ 
sented in Section 10.6, we must compute eleven matrix integrals which can 
be expressed as one of the following two forms: 


G.(t) - / e A(t s) Be Cs ds (10-46) 

1 0 

or 

H(t) - / / e A(t_8) Be C(s " r) De Er dr ds (10-47) 

0 0 


where the (constant) matrices A, B, C, D, and E are specified. This 
problem, and its solution, are analogous to that discussed in the 
preceding subsection. In order to compute Gj and Hj, we define the 
following matrix exponential: 

T(t) = r(A, B, C, D, E, t) - e Vt (10-48) 

where v is the constant matrix [Ref. (10-40)] 


v 


A 


0 

L o 


B 0-1 } nj 

C tt )n 2 

0 E J } n 3 

"2 "3 
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As shown in Reference 10-40, T has the representation: 


Gj(t) 

HjCt) 

F 2 (t) 

G 2 (t) 

0 

F 3 (t) 

»2 

"3 


which includes the desired integrals Gj and In addition, 


Using the semi-group properties of matrix exponentials, the partitions of 
r can be recursively generated from the following formulas: 

F 1 (t + At) - F i (At)F i (t) , F 1 (0) - I , i - 1, 2, 3 

G A (t + At) = F 1 (At)G 1 (t) + G 1 (At)F 1+1 (t), G^O) - 0 , i - 1, 2 

H^t + At) - FjCAOHjCt) + G 1 (At)G 2 (t) + H 1 (At)F 3 (t) , H^O) - 0 

When the recursion relationships above are post-multiplied by vectors, 
the following modified propagation equations result: 


IVV.\V* V.'iA'V-V'> - -■ 













x^t + At) = F^t + At)x Q = F^(At)x^(t) , i * 1, 2, 3 

(10-49) 

x ± (t + At) - Gj_^(t + At)x Q = F 13 (At)x i (t) + G 13 (At)x i 2 (t) , i - 4, 5 

(10-50) 


At + At) = H,(t + At)x. = F (At)x.(t) + G (At)x (t) + H.(At)x_(t) 


where 



i = 1 , 2, 3 
i = 4, 5, 6 


The vector recursion relationships above are particularly useful 
for generating the integral solution for x r (t) in Eq. (10-24) at 
discrete times, since matrix-vector multiplications are much less costly 
than matrix-matrix multiplications. Furthermore, in generating 
matrix-vector products involving matrices of the form in Eq. (10-46), 
only X 2 and X 4 in Eqs. (10-49), (10-50) must be propagated. 

Similarly, in generating matrix-vector products involving matrices of the 
form in Eq. (10-47), only X3, X5 and xg must be propagated. 

10.10 Conclusions 

Closed-form solutions for the feedback gains, state trajectory, 
and residual mode response have been developed for the tracking/ 
disturbance-accommodation controller. This has been made possible by the 
fact that the Fourier series representation of the disturbances may be 
expressed in an exponential matrix form. In addition, recursion rela¬ 
tionships have been developed for propagating the gains, state trajectory 
and residual mode response in a computationally efficient way. Numerical 
results of several test case maneuvers demonstrate the utility of the 




254 







formulations and show that the tracking/disturbance-accommodation con¬ 
troller outperforms the simple linear tracking controller which does not 
use Information about the disturbances. The test case examples also show 
the effect that control-rate smoothing has on smoothing out the control 
time history for the linear tracker. 

A computationally efficient algorithm has also been presented for 
obtaining the complex Fourier transform of a class of vector functions 
which frequently occur in modern control theory. The basic algorithm 
requires: (1) evaluation of a single matrix exponential for the dynamics 
of the time-varying control; (2) solution for either the right and left 
eigenvectors or a real Schur decomposition of the constant control 
dynamics matrix; (3) sequential solution for the reducing subspace 
transformation vector; and (4) evaluation of a single scalar complex 
exponential. A low order example application is shown to demonstrate the 
utility of the above formulation. 
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